Uncertainty Quantification Group

Available Links

• https://arxiv.org/abs/2411.12726

Abstract

We present LazyDINO, a transport map variational inference method for fast, scalable, and efficiently amortized solutions of high-dimensional nonlinear Bayesian inverse problems with expensive parameter-to-observable (PtO) maps. Our method consists of an offline phase in which we construct a derivative-informed neural surrogate of the PtO map using joint samples of the PtO map and its Jacobian. During the online phase, when given observational data, we seek rapid posterior approximation using surrogate-driven training of a lazy map [Brennan et al., NeurIPS, (2020)], i.e., a structure-exploiting transport map with low-dimensional nonlinearity. The trained lazy map then produces approximate posterior samples or density evaluations. Our surrogate construction is optimized for amortized Bayesian inversion using lazy map variational inference. We show that (i) the derivative-based reduced basis architecture [O'Leary-Roseberry et al., Comput. Methods Appl. Mech. Eng., 388 (2022)] minimizes the upper bound on the expected error in surrogate posterior approximation, and (ii) the derivative-informed training formulation [O'Leary-Roseberry et al., J. Comput. Phys., 496 (2024)] minimizes the expected error due to surrogate-driven transport map optimization. Our numerical results demonstrate that LazyDINO is highly efficient in cost amortization for Bayesian inversion. We observe one to two orders of magnitude reduction of offline cost for accurate posterior approximation, compared to simulation-based amortized inference via conditional transport and conventional surrogate-driven transport. In particular, LazyDINO outperforms Laplace approximation consistently using fewer than 1000 offline samples, while other amortized inference methods struggle and sometimes fail at 16,000 offline samples.

BibTeX Citation

@article{ cao-lazydino-2026,
  keywords = { transport,dimension reduction,surrogates,amortized inference },
  title = { LazyDINO: Fast, scalable, and efficiently amortized Bayesian inversion via structure-exploiting and surrogate-driven measure transport },
  author = { L. Cao and J. Chen and M. Brennan and T. O'Leary-Roseberry and Y. M. Marzouk and O. Ghattas },
  journal = { Journal of Machine Learning Research },
  volume = { 27 },
  number = { 24 },
  pages = { 1--71 },
  month = { 2 },
  year = { 2026 },
  editor = {  },
}

Available Links

• https://arxiv.org/abs/2411.08039

Abstract

BibTeX Citation

@article{ lu-fluidleakagefaultinstability-2025,
  title = { Uncertainty Quantification of Fluid Leakage and Fault Instability in Geologic CO₂ Storage },
  author = { H. Lu and L. Saló--Salgado and Y. M. Marzouk and R. Juanes },
  journal = { Water Resources Research },
  volume = { 61 },
  number = { 10 },
  year = { 2025 },
  publisher = { Wiley Online Library },
  editor = {  },
}

Available Links

• https://arxiv.org/abs/2307.12438

Abstract

We introduce a multifidelity estimator of covariance matrices formulated as the solution to a regression problem on the manifold of symmetric positive definite matrices. The estimator is positive definite by construction, and the Mahalanobis distance minimized to obtain it possesses properties which enable practical computation. We show that our manifold regression multifidelity (MRMF) covariance estimator is a maximum likelihood estimator under a certain error model on manifold tangent space. More broadly, we show that our Riemannian regression framework encompasses existing multifidelity covariance estimators constructed from control variates. We demonstrate via numerical examples that our estimator can provide significant decreases, up to one order of magnitude, in squared estimation error relative to both single-fidelity and other multifidelity covariance estimators. Furthermore, preservation of positive definiteness ensures that our estimator is compatible with downstream tasks, such as data assimilation and metric learning, in which this property is essential.

BibTeX Citation

@article{ maurais-covarianceregression-2025,
  keywords = { covariance, regression, SPD },
  doi = { doi:10.1137/23M159247X },
  title = { Multifidelity Covariance Estimation via Regression on the Manifold of Symmetric Positive Definite Matrices },
  author = { A. Maurais and T. Alsup and B. Peherstorfer and Y. M. Marzouk },
  journal = { SIAM Journal on Mathematics of Data Science },
  volume = { 7 },
  number = { 1 },
  pages = { 189-223 },
  month = { 1 },
  year = { 2025 },
  editor = {  },
}

Available Links

• https://projecteuclid.org/journals/bernoulli/volume-30/issue-4/Principal-feature-detection-via-%CF%95-Sobolev-inequalities/10.3150/23-BEJ1702.full

Abstract

We investigate the approximation of high-dimensional target measures as low-dimensional updates of a dominating reference measure. This approximation class replaces the associated density with the composition of: (i) a feature map that identifies the leading principal components or features of the target measure, relative to the reference, and (ii) a low-dimensional profile function. When the reference measure satisfies a subspace ϕ-Sobolev inequality, we construct a computationally tractable approximation that yields certifiable error guarantees with respect to the Amari α-divergences. Our construction proceeds in two stages. First, for any feature map and any α-divergence, we obtain an analytical expression for the optimal profile function. Second, for linear feature maps, the principal features are obtained from eigenvectors of a matrix involving gradients of the log-density. Neither step requires explicit access to normalizing constants. Notably, by leveraging the ϕ-Sobolev inequalities, we demonstrate that these features universally certify approximation errors across the range of α-divergences α ∈(0,1]. We then propose an application to Bayesian inverse problems and provide an analogous construction with approximation guarantees that hold in expectation over the data. We conclude with an extension of the proposed dimension reduction strategy to nonlinear feature maps.

BibTeX Citation

@article{ li-featuredetection-2024,
  keywords = { Amari α-divergences; Bayesian inference; feature detection; gradient-based dimension reduction; principal components; φ-Sobolev inequalities },
  doi = { 10.3150/23-BEJ1702 },
  title = { Principal feature detection via ϕ-Sobolev inequalities },
  author = { M. T. C. Li and Y. M. Marzouk and O. Zahm },
  journal = { Bernoulli },
  volume = { 30 },
  number = { 4 },
  pages = { 2979-3003 },
  month = { 0 },
  year = { 2024 },
  editor = {  },
}

Available Links

• https://www.cambridge.org/core/journals/acta-numerica/article/optimal-experimental-design-formulations-and-computations/38BBD0DC1A0386FDF306B6C0167DF7D9

Abstract

Questions of ‘how best to acquire data’ are essential to modelling and prediction in the natural and social sciences, engineering applications, and beyond. Optimal experimental design (OED) formalizes these questions and creates computational methods to answer them. This article presents a systematic survey of modern OED, from its foundations in classical design theory to current research involving OED for complex models. We begin by reviewing criteria used to formulate an OED problem and thus to encode the goal of performing an experiment. We emphasize the flexibility of the Bayesian and decision-theoretic approach, which encompasses information-based criteria that are well-suited to nonlinear and non-Gaussian statistical models. We then discuss methods for estimating or bounding the values of these design criteria; this endeavour can be quite challenging due to strong nonlinearities, high parameter dimension, large per-sample costs, or settings where the model is implicit. A complementary set of computational issues involves optimization methods used to find a design; we discuss such methods in the discrete (combinatorial) setting of observation selection and in settings where an exact design can be continuously parametrized. Finally we present emerging methods for sequential OED that build non-myopic design policies, rather than explicit designs; these methods naturally adapt to the outcomes of past experiments in proposing new experiments, while seeking coordination among all experiments to be performed. Throughout, we highlight important open questions and challenges.

BibTeX Citation

@article{ huan-oedreview-2024,
  keywords = { optimal experimental design },
  doi = { 10.1017/S0962492924000023 },
  title = { Optimal experimental design: Formulations and computations },
  author = { X. Huan and J. Jagalur and Y. M. Marzouk },
  journal = { Acta Numerica },
  volume = { 33 },
  pages = { 715-840 },
  month = { 7 },
  year = { 2024 },
  editor = {  },
}

Available Links

• https://arxiv.org/pdf/2206.11343

Abstract

We consider the Bayesian calibration of models describing the phenomenon of block copolymer (BCP) self-assembly, which is to infer model parameters and their uncertainty from noisy image data produced by microscopy or X-ray scattering characterization of BCP structures. To account for the long-range disorder in BCP structures, we introduce auxiliary variables in the model to represent this aleatory uncertainty. These variables, however, result in an integrated likelihood function for high-dimensional image data that is generally intractable to evaluate. We tackle this Bayesian inference problem using a likelihood-free inference (LFI) approach based on measure transport together with the construction of summary statistics for the image data. We show that expected information gains (EIG) from observed data about the model parameters can be computed with no significant additional cost. Lastly, we present a numerical case study using the Ohta–Kawasaki model for diblock copolymer thin film self-assembly and top-down microscopy characterization. We introduce several domain-specific energy and Fourier-based summary statistics for calibration and quantify their informativeness using EIG. The effect of various data corruptions, summary statistics, and experimental designs on the calibration results are studied using the proposed LFI method.

BibTeX Citation

@article{ baptista-bayesiancopolymer-2024,
  keywords = { Block copolymers, Uncertainty quantification, Likelihood-free inference, Measure transport, Expected information gain, Ohta–Kawasaki model },
  doi = { 10.1016/j.jcp.2024.112844 },
  title = { Bayesian model calibration for block copolymer self-assembly: Likelihood-free inference and expected information gain computation via measure transport },
  author = { R. Baptista and L. Cao and J. Chen and O. Ghattas and F. Y. Li and Y. M. Marzouk and J. T. Oden },
  journal = { Journal of Computational Physics },
  volume = { 503 },
  number = { 112844 },
  month = { 4 },
  year = { 2024 },
  editor = {  },
}

Available Links

• https://www.jmlr.org/papers/volume25/23-1280/23-1280.pdf

Abstract

Ordinary differential equations (ODEs), via their induced flow maps, provide a powerful framework to parameterize invertible transformations for representing complex probability distributions. While such models have achieved enormous success in machine learning, little is known about their statistical properties. This work establishes the first general nonparametric statistical convergence analysis for distribution learning via ODE models trained through likelihood maximization. We first prove a convergence theorem applicable to arbitrary velocity field classes F satisfying certain simple boundary constraints. This general result captures the trade-off between the approximation error and complexity of the ODE model. We show that the latter can be quantified via the C1-metric entropy of the class F. We then apply this general framework to the setting of Ck-smooth target densities, and establish nearly minimax-optimal convergence rates for two relevant velocity field classes F: Ck functions and neural networks. The latter is the practically important case of neural ODEs. Our results also provide insight on how the choice of velocity field class, and the dependence of this choice on sample size (e.g., the scaling of neural network classes), impact statistical performance.

BibTeX Citation

@article{ marzouk-distributionneuralode-2024,
  keywords = { Neural differential equations, normalizing flows, density estimation, nonparametric statistics, M-estimation },
  title = { Distribution Learning via Neural Differential Equations: A Nonparametric Statistical Perspective },
  author = { Y. M. Marzouk and Z. R. Ren and S. Wang and J. Zech },
  journal = { Journal of Machine Learning Research },
  volume = { 25 },
  number = { 232 },
  pages = { 1-61 },
  month = { 6 },
  year = { 2024 },
  editor = {  },
}

Available Links

• https://www.jmlr.org/papers/volume25/21-0022/21-0022.pdf

Abstract

Undirected probabilistic graphical models represent the conditional dependencies, or Markov properties, of a collection of random variables. Knowing the sparsity of such a graphical model is valuable for modeling multivariate distributions and for efficiently performing inference. While the problem of learning graph structure from data has been studied extensively for certain parametric families of distributions, most existing methods fail to consistently recover the graph structure for non-Gaussian data. Here we propose an algorithm for learning the Markov structure of continuous and non-Gaussian distributions. To characterize conditional independence, we introduce a score based on integrated Hessian information from the joint log-density, and we prove that this score upper bounds the conditional mutual information for a general class of distributions. To compute the score, our algorithm SING estimates the density using a deterministic coupling, induced by a triangular transport map, and iteratively exploits sparse structure in the map to reveal sparsity in the graph. For certain non-Gaussian datasets, we show that our algorithm recovers the graph structure even with a biased approximation to the density. Among other examples, we apply SING to learn the dependencies between the states of a chaotic dynamical system with local interactions.

BibTeX Citation

@article{ baptista-nongaussiangraphical-2024,
  keywords = { Undirected graphical models, structure learning, non-Gaussian distributions, conditional mutual information, transport map, sparsity },
  title = { Learning Non-Gaussian Graphical Models via Hessian Scores and Triangular Transport },
  author = { R. Baptista and R. Morrison and O. Zahm and Y. M. Marzouk },
  journal = { Journal of Machine Learning Research },
  volume = { 25 },
  number = { 85 },
  pages = { 1-46 },
  month = { 0 },
  year = { 2024 },
  editor = {  },
}

Available Links

• https://arxiv.org/pdf/2305.03103

Abstract

Optimization is a key tool for scientific and engineering applications; however, in the presence of models affected by uncertainty, the optimization formulation needs to be extended to consider statistics of the quantity of interest. Optimization under uncertainty (OUU) deals with this endeavor and requires uncertainty quantification analyses at several design locations; i.e., its overall computational cost is proportional to the cost of performing a forward uncertainty analysis at each design location. An OUU workflow has two main components: an inner loop strategy for the computation of statistics of the quantity of interest, and an outer loop optimization strategy tasked with finding the optimal design, given a merit function based on the inner loop statistics. In this work, we propose to alleviate the cost of the inner loop uncertainty analysis by leveraging the so-called multilevel Monte Carlo (MLMC) method, which is able to allocate resources over multiple models with varying accuracy and cost. The resource allocation problem in MLMC is formulated by minimizing the computational cost given a target variance for the estimator. We consider MLMC estimators for statistics usually employed in OUU workflows and solve the corresponding allocation problem. For the outer loop, we consider a derivative-free optimization strategy implemented in the SNOWPAC library; our novel strategy is implemented and released in the Dakota software toolkit. We discuss several numerical test cases to showcase the features and performance of our approach with respect to its Monte Carlo single fidelity counterpart.

BibTeX Citation

@article{ menhorn-multileveloptimization-2024,
  keywords = { optimization under uncertainty, multilevel Monte Carlo, uncertainty quantification },
  doi = { 10.1615/Int.J.UncertaintyQuantification.2023048049 },
  title = { Multilevel Monte Carlo estimators for derivative-free optimization under uncertainty },
  author = { F. Menhorn and G. Geraci and D. T. Seidl and Y. M. Marzouk and M. S. Eldred and H.-J. Bungartz },
  journal = { International Journal for Uncertainty Quantification },
  volume = { 14 },
  number = { 3 },
  pages = { 21-65 },
  month = { 11 },
  year = { 2023 },
  editor = {  },
}

Available Links

• https://arxiv.org/pdf/2310.16975

Abstract

We present two neural network approaches that approximate the solutions of static and dynamic conditional optimal transport (COT) problems. Both approaches enable conditional sampling and conditional density estimation, which are core tasks in Bayesian inference–particularly in the simulation-based ("likelihood-free") setting. Our methods represent the target conditional distributions as transformations of a tractable reference distribution and, therefore, fall into the framework of measure transport. Although many measure transport approaches model the transformation as COT maps, obtaining the map is computationally challenging, even in moderate dimensions. To improve scalability, our numerical algorithms use neural networks to parameterize COT maps and further exploit the structure of the COT problem. Our static approach approximates the map as the gradient of a partially input-convex neural network. It uses a novel numerical implementation to increase computational efficiency compared to state-of-the-art alternatives. Our dynamic approach approximates the conditional optimal transport via the flow map of a regularized neural ODE; compared to the static approach, it is slower to train but offers more modeling choices and can lead to faster sampling. We demonstrate both algorithms numerically, comparing them with competing state-of-the-art approaches, using benchmark datasets and simulation-based Bayesian inverse problems.

BibTeX Citation

@article{ wang-conditionalot-2025,
  keywords = { measure transport, generative modeling, optimal transport, Bayesian inference, inverse problems, uncertainty quantification },
  doi = { 10.1137/24M1678659 },
  title = { Efficient Neural Network Approaches for Conditional Optimal Transport with Applications in Bayesian Inference },
  author = { Z. O. Wang and R. Baptista and Y. M. Marzouk and L. Ruthotto and D. Verma },
  journal = { SIAM Journal on Scientific Computing },
  volume = { 47 },
  number = { 4 },
  pages = { C979--C1005 },
  month = { 8 },
  year = { 2025 },
  editor = {  },
}

Available Links

• https://academic.oup.com/gji/article/234/2/1169/7081312

Abstract

A velocity model is generally an imperfect representation of the subsurface, which cannot precisely account for the 3-D inhomogeneities of Earth structure. We present a Bayesian moment tensor inversion framework for applications where reliable, tomography-based, velocity model reconstructions are not available. In particular, synthetic data generated using a 3-D model (SEG-EAGE Overthrust) are inverted using a layered medium model. We use a likelihood function derived from an optimal transport distance—specifically, the transport-Lagrangian distance introduced by Thorpe et al.—and show that this formulation yields inferences that are robust to misspecification of the velocity model. We establish several quantitative metrics to evaluate the performance of the proposed Bayesian framework, comparing it to Bayesian inversion with a standard Gaussian likelihood. We also show that the non-double-couple component of the recovered mechanisms drastically diminishes when the impact of velocity model misspecification is mitigated.

BibTeX Citation

@article{ scarinci-tensorinversion-2023,
  keywords = { Inverse theory; Probability distributions; Statistical methods; Earthquake source observations; Induced seismicity; Waveform inversion },
  doi = { 10.1093/gji/ggad116 },
  title = { Robust Bayesian moment tensor inversion with optimal transport misfits: layered medium approximations to the 3-D SEG-EAGE overthrust velocity model },
  author = { A. Scarinci and U. bin Waheed and C. Gu and X. Ren and B. M. Dia and S. Kaka and M. Fehler and Y. M. Marzouk },
  journal = { Geophysical Journal International },
  volume = { 234 },
  number = { 2 },
  pages = { 1169-1190 },
  month = { 8 },
  year = { 2023 },
  editor = {  },
}

Available Links

• https://arxiv.org/abs/2210.17000

Abstract

Smoothers are algorithms for Bayesian time series re-analysis. Most operational smoothers rely either on affine Kalman-type transformations or on sequential importance sampling. These strategies occupy opposite ends of a spectrum that trades computational efficiency and scalability for statistical generality and consistency: non-Gaussianity renders affine Kalman updates inconsistent with the true Bayesian solution, while the ensemble size required for successful importance sampling can be prohibitive. This paper revisits the smoothing problem from the perspective of measure transport, which offers the prospect of consistent prior-to-posterior transformations for Bayesian inference. We leverage this capacity by proposing a general ensemble framework for transport-based smoothing. Within this framework, we derive a comprehensive set of smoothing recursions based on nonlinear transport maps and detail how they exploit the structure of state-space models in fully non-Gaussian settings. We also describe how many standard Kalman-type smoothing algorithms emerge as special cases of our framework. A companion paper explores the implementation of nonlinear ensemble transport smoothers in greater depth.

BibTeX Citation

@article{ ramgraber-smoothingpartone-2022,
  keywords = { Data assimilation, smoothing, ensemble methods, triangular transport },
  doi = { 10.1016/j.jcpx.2023.100134 },
  title = { Ensemble transport smoothing. Part 1: unified framework },
  author = { M. Ramgraber and R. Baptista and D. McLaughlin and Y. M. Marzouk },
  journal = { Journal of Computational Physics: X },
  volume = { 17 },
  number = { 100134 },
  month = { 11 },
  year = { 2023 },
  editor = {  },
}

Available Links

• https://arxiv.org/abs/2112.00713

Abstract

Bayesian inference provides a systematic means of quantifying uncertainty in the solution of the inverse problem. However, solution of Bayesian inverse problems governed by complex forward models described by partial differential equations (PDEs) remains prohibitive with black-box Markov chain Monte Carlo (MCMC) methods. We present hIPPYlib-MUQ, an extensible and scalable software framework that contains implementations of state-of-the art algorithms aimed to overcome the challenges of high-dimensional, PDE-constrained Bayesian inverse problems. hIPPYlib-MUQ integrates two complementary open-source software packages. hIPPYlib solves PDE-constrained inverse problems using automatically-generated adjoint-based derivatives, but it lacks full Bayesian capabilities. MUQ provides numerous powerful Bayesian inversion algorithms, but expects forward models to come equipped with derivatives to permit large-scale solution. By combining these two libraries, we created a robust, scalable, and efficient software framework that can be used to tackle complex large-scale Bayesian inverse problems across a broad spectrum of scientific and engineering disciplines. To illustrate the capabilities of hIPPYlib-MUQ, we compare a number of MCMC methods on several high-dimensional Bayesian inverse problems. The results demonstrate that large (∼50×) speedups over conventional black box and gradient-based MCMC algorithms can be obtained by exploiting Hessian information (from the log-posterior), underscoring the power of the integrated hIPPYlib-MUQ framework.

BibTeX Citation

@article{ kim-hippylibmuq-2021,
  doi = { 10.48550/arXiv.2112.00713 },
  title = { hIPPYlib-MUQ: a Bayesian inference software framework for integration of data with complex predictive models under uncertainty },
  author = { K.-T. Kim and U. Villa and M. Parno and Y. M. Marzouk and O. Ghattas and N. Petra },
  journal = { ACM Transactions on Mathematical Software },
  month = { 0 },
  year = { 2022 },
  editor = {  },
}

Available Links

• https://arxiv.org/abs/1907.00389
• https://doi.org/10.1137/20M1312204

Abstract

We consider filtering in high-dimensional non-Gaussian state-space models with intractable transition kernels, nonlinear and possibly chaotic dynamics, and sparse observations in space and time. We propose a novel filtering methodology that harnesses transportation of measures, convex optimization, and ideas from probabilistic graphical models to yield robust ensemble approximations of the filtering distribution in high dimensions. Our approach can be understood as the natural generalization of the ensemble Kalman filter (EnKF) to nonlinear updates, using stochastic or deterministic couplings. The use of nonlinear updates can reduce the intrinsic bias of the EnKF at a marginal increase in computational cost. We avoid any form of importance sampling and introduce non-Gaussian localization approaches for dimension scalability. Our framework achieves state-of-the-art tracking performance on challenging configurations of the Lorenz-96 model in the chaotic regime.

BibTeX Citation

@article{ art_79,
  keywords = { nonlinear filtering, state-space models, couplings, transport maps, ensemble Kalman filter, graphical models, localization, approximate Bayesian computation },
  doi = { 10.1137/20M1312204 },
  title = { Coupling techniques for nonlinear ensemble filtering },
  author = { A. Spantini and R. Baptista and Y. M. Marzouk },
  journal = { SIAM Review },
  volume = { 64 },
  number = { 4 },
  pages = { 921–953 },
  month = { 0 },
  year = { 2022 },
  editor = {  },
}

Available Links

• https://arxiv.org/abs/2111.07164
• https://doi.org/10.1002/nla.2467

Abstract

Very often, in the course of uncertainty quantification tasks or data analysis, one has to deal with high-dimensional random variables. Here the interest is mainly to compute characterizations like the entropy, the Kullback–Leibler divergence, more general -divergences, or other such characteristics based on the probability density. The density is often not available directly, and it is a computational challenge to just represent it in a numerically feasible fashion in case the dimension is even moderately large. It is an even stronger numerical challenge to then actually compute said characteristics in the high-dimensional case. In this regard it is proposed to approximate the discretized density in a compressed form, in particular by a low-rank tensor. This can alternatively be obtained from the corresponding probability characteristic function, or more general representations of the underlying random variable. The mentioned characterizations need point-wise functions like the logarithm. This normally rather trivial task becomes computationally difficult when the density is approximated in a compressed resp. low-rank tensor format, as the point values are not directly accessible. The computations become possible by considering the compressed data as an element of an associative, commutative algebra with an inner product, and using matrix algorithms to accomplish the mentioned tasks. The representation as a low-rank element of a high order tensor space allows to reduce the computational complexity and storage cost from exponential in the dimension to almost linear.

BibTeX Citation

@article{ litvinenko-tensordiv-2022,
  doi = { 10.1002/nla.2467 },
  title = { Computing f-divergences and distances of high-dimensional probability density functions },
  author = { A. Litvinenko and Y. M. Marzouk and H. G. Matthies and M. Scavino and A. Spantini },
  journal = { Numerical Linear Algebra with Applications },
  volume = { in press },
  month = { 0 },
  year = { 2022 },
  editor = {  },
}

Available Links

• https://arxiv.org/abs/2108.08247
• https://doi.org/10.1007/s11222-022-10147-6

Abstract

We introduce a novel geometry-informed irreversible perturbation that accelerates convergence of the Langevin algorithm for Bayesian computation. It is well documented that there exist perturbations to the Langevin dynamics that preserve its invariant measure while accelerating its convergence. Irreversible perturbations and reversible perturbations (such as Riemannian manifold Langevin dynamics (RMLD)) have separately been shown to improve the performance of Langevin samplers. We consider these two perturbations simultaneously by presenting a novel form of irreversible perturbation for RMLD that is informed by the underlying geometry. Through numerical examples, we show that this new irreversible perturbation can improve estimation performance over irreversible perturbations that do not take the geometry into account. Moreover we demonstrate that irreversible perturbations generally can be implemented in conjunction with the stochastic gradient version of the Langevin algorithm. Lastly, while continuous-time irreversible perturbations cannot impair the performance of a Langevin estimator, the situation can sometimes be more complicated when discretization is considered. To this end, we describe a discrete-time example in which irreversibility increases both the bias and variance of the resulting estimator.

BibTeX Citation

@article{ zhang-girr-2021,
  doi = { 10.1007/s11222-022-10147-6 },
  title = { Geometry-informed irreversible perturbations for accelerated convergence of Langevin dynamics },
  author = { B. J. Zhang and Y. M. Marzouk and K. Spiliopoulos },
  journal = { Statistics and Computing },
  volume = { 32 },
  number = { 78 },
  month = { 0 },
  year = { 2022 },
  editor = {  },
}

Available Links

• https://arxiv.org/abs/2107.13422
• https://doi.org/10.1007/s00365-022-09570-9

Abstract

For two probability measures $\rho$ and $\pi$ on $[−1,1]^\mathbb{N}$ we investigate the approximation of the triangular Knothe-Rosenblatt transport $T: [−1,1]^\mathbb{N} \to [−1,1]^\mathbb{N}$ that pushes forward $\rho$ to $\pi$. Under suitable assumptions, we show that $T$ can be approximated by rational functions without suffering from the curse of dimension. Our results are applicable to posterior measures arising in certain inference problems where the unknown belongs to an (infinite-dimensional) Banach space. In particular, we show that it is possible to efficiently approximately sample from certain high-dimensional measures by transforming a lower-dimensional latent variable.

BibTeX Citation

@article{ zech-infdim-2021,
  doi = { 10.1007/s00365-022-09570-9 },
  title = { Sparse approximation of triangular transports. Part II: the infinite dimensional case },
  author = { J. Zech and Y. M. Marzouk },
  journal = { Constructive Approximation },
  volume = { 55 },
  pages = { 987–1036 },
  month = { 0 },
  year = { 2022 },
  editor = {  },
}

Available Links

• https://arxiv.org/abs/2101.07330
• https://doi.org/10.1016/j.jcp.2022.111025

Abstract

We exploit the relationship between the stochastic Koopman operator and the Kolmogorov backward equation to construct importance sampling schemes for stochastic differential equations. Specifically, we propose using eigenfunctions of the stochastic Koopman operator to approximate the Doob transform for an observable of interest (e.g., associated with a rare event) which in turn yields an approximation of the corresponding zero-variance importance sampling estimator. Our approach is broadly applicable and systematic, treating non-normal systems, non-gradient systems, and systems with oscillatory dynamics or rank-deficient noise in a common framework. In nonlinear settings where the stochastic Koopman eigenfunctions cannot be derived analytically, we use dynamic mode decomposition (DMD) methods to approximate them numerically, but the framework is agnostic to the particular numerical method employed. Numerical experiments demonstrate that even coarse approximations of a few eigenfunctions, where the latter are built from non-rare trajectories, can produce effective importance sampling schemes for rare events.

BibTeX Citation

@article{ zhang-kfr-2021,
  doi = { 10.1016/j.jcp.2022.111025 },
  title = { A Koopman framework for rare event simulation in stochastic differential equations },
  author = { B. J. Zhang and T. Sahai and Y. M. Marzouk },
  journal = { Journal of Computational Physics },
  volume = { 456 },
  pages = { 111025 },
  month = { 0 },
  year = { 2022 },
  editor = {  },
}

Available Links

• https://arxiv.org/abs/2004.12102

Abstract

We propose a differential geometric approach for building families of low-rank covariance matrices, via interpolation on low-rank matrix manifolds. In contrast with standard parametric covariance classes, these families offer significant flexibility for problem-specific tailoring via the choice of “anchor” matrices for interpolation, for instance over a grid of relevant conditions describing the underlying stochastic process. The interpolation is computationally tractable in high dimensions, as it only involves manipulations of low-rank matrix factors. We also consider the problem of covariance identification, i.e., selecting the most representative member of the covariance family given a data set. In this setting, standard procedures such as maximum likelihood estimation are nontrivial because the covariance family is rank-deficient; we resolve this issue by casting the identification problem as distance minimization. We demonstrate the utility of these differential geometric families for interpolation and identification in a practical application: wind field covariance approximation for unmanned aerial vehicle navigation.

BibTeX Citation

@article{ musolas-lowrank-2020,
  doi = { https://doi.org/10.1007/s10543-021-00867-y },
  title = { Low-rank multi-parametric covariance identification },
  author = { A. Musolas and E. Massart and J. Hendrickx and P.-A. Absil and Y. M. Marzouk },
  journal = { BIT Numerical Mathematics },
  volume = { 62 },
  pages = { 221–249 },
  month = { 0 },
  year = { 2022 },
  editor = {  },
}

Available Links

• https://gmd.copernicus.org/articles/14/4319/2021/

Abstract

Estimating parameters of chaotic geophysical models is challenging due to their inherent unpredictability. These models cannot be calibrated with standard least squares or filtering methods if observations are temporally sparse. Obvious remedies, such as averaging over temporal and spatial data to characterize the mean behavior, do not capture the subtleties of the underlying dynamics. We perform Bayesian inference of parameters in high-dimensional and computationally demanding chaotic dynamical systems by combining two approaches: (i) measuring model–data mismatch by comparing chaotic attractors and (ii) mitigating the computational cost of inference by using surrogate models. Specifically, we construct a likelihood function suited to chaotic models by evaluating a distribution over distances between points in the phase space; this distribution defines a summary statistic that depends on the geometry of the attractor, rather than on pointwise matching of trajectories. This statistic is computationally expensive to simulate, compounding the usual challenges of Bayesian computation with physical models. Thus, we develop an inexpensive surrogate for the log likelihood with the local approximation Markov chain Monte Carlo method, which in our simulations reduces the time required for accurate inference by orders of magnitude. We investigate the behavior of the resulting algorithm with two smaller-scale problems and then use a quasi-geostrophic model to demonstrate its large-scale application.

BibTeX Citation

@article{ springer-ebi-2021,
  doi = { https://doi.org/10.5194/gmd-14-4319-2021 },
  title = { Efficient Bayesian inference for large chaotic dynamical systems },
  author = { S. Springer and H. Haario and J. Susiluoto and A. Bibov and A. Davis and Y. M. Marzouk },
  journal = { Geoscientific Model Development },
  volume = { 14 },
  number = { 7 },
  pages = { 4319–4333 },
  month = { 0 },
  year = { 2021 },
  editor = {  },
}

Available Links

• https://arxiv.org/abs/2001.01805

Abstract

Statistical modeling of spatiotemporal phenomena often requires selecting a covariance matrix from a covariance class. Yet standard parametric covariance families can be insufficiently flexible for practical applications, while nonparametric approaches may not easily allow certain kinds of prior knowledge to be incorporated. We propose instead to build covariance families out of geodesic curves. These covariances offer more flexibility for problem-specific tailoring than classical parametric families and are preferable to simple convex combinations. Once the covariance family has been chosen, one typically needs to select a representative member by solving an optimization problem, e.g., by maximizing the likelihood of a data set. We consider instead a differential geometric interpretation of this problem: minimizing the geodesic distance to a sample covariance matrix (``natural projection''). Our approach is consistent with the notion of distance employed to build the covariance family and does not require assuming a particular probability distribution for the data. We show that natural projection and maximum likelihood estimation within the covariance family are locally equivalent up to second order. We also demonstrate that natural projection may yield more accurate estimates with noise-corrupted data.

BibTeX Citation

@article{ musolas-fullrank-2020,
  doi = { https://doi.org/10.1137/19M1284646 },
  title = { Geodesically parameterized covariance estimation },
  author = { A. Musolas and S. T. Smith and Y. M. Marzouk },
  journal = { SIAM Journal on Matrix Analysis and Applications },
  volume = { 42 },
  number = { 2 },
  pages = { 528–556 },
  month = { 0 },
  year = { 2021 },
  editor = {  },
}

Available Links

• https://arxiv.org/abs/2003.07991
• https://doi.org/10.1088/1361-6420/abb2fa

Abstract

This paper suggests a framework for the learning of discretizations of expensive forward models in Bayesian inverse problems. The main idea is to incorporate the parameters governing the discretization as part of the unknown to be estimated within the Bayesian machinery. We numerically show that in a variety of inverse problems arising in mechanical engineering, signal processing and the geosciences, the observations contain useful information to guide the choice of discretization.

BibTeX Citation

@article{ bigoni-discretization-2020,
  doi = { 10.1088/1361-6420/abb2fa },
  title = { Data-driven forward discretizations for Bayesian inversion },
  author = { D. Bigoni and Y. Chen and N. Garcia-Trillos and Y. M. Marzouk and D. Sanz-Alonso },
  journal = { Inverse Problems },
  volume = { 36 },
  number = { 10 },
  pages = { 105008 },
  month = { 0 },
  year = { 2020 },
  editor = {  },
}

Available Links

• https://arxiv.org/abs/2006.06798
• https://doi.org/10.1088/1741-4326/abae85

Abstract

We present a fully Bayesian approach for the inference of radial profiles of impurity transport coefficients and compare its results to neoclassical, gyrofluid and gyrokinetic modeling. Using nested sampling, the Bayesian Impurity Transport InferencE (BITE) framework can handle complex parameter spaces with multiple possible solutions, offering great advantages in interpretative power and reliability with respect to previously demonstrated methods. BITE employs a forward model based on the pySTRAHL package, built on the success of the well-known STRAHL code [Dux, IPP Report, 2004], to simulate impurity transport in magnetically-confined plasmas. In this paper, we focus on calcium (Ca, Z=20) Laser Blow-Off injections into Alcator C-Mod plasmas. Multiple Ca atomic lines are diagnosed via high-resolution X-ray Imaging Crystal Spectroscopy and Vacuum Ultra-Violet measurements. We analyze a sawtoothing I-mode discharge for which neoclassical and turbulent (quasilinear and nonlinear) predictions are also obtained. We find good agreement in diffusion across the entire radial extent, while turbulent convection and density profile peaking are estimated to be larger in experiment than suggested by theory. Efforts and challenges associated with the inference of experimental pedestal impurity transport are discussed.

BibTeX Citation

@article{ sciortino-impurity-2020,
  doi = { 10.1088/1741-4326/abae85 },
  title = { Inference of experimental radial impurity transport on Alcator C-Mod: Bayesian parameter estimation and model selection },
  author = { F. Sciortino and N. T. Howard and E. S. Marmar and T. Odstrcil and N. M. Cao and R. Dux and A. E. Hubbard and J. W. Hughes and J. H. Irby and Y. M. Marzouk and L. M. Milanese },
  journal = { Nuclear Fusion },
  volume = { 60 },
  pages = { 126014 },
  month = { 0 },
  year = { 2020 },
  editor = {  },
}

Available Links

• https://arxiv.org/abs/1906.10098

Abstract

Acoustic emission (AE) is a widely used technology to study source mechanisms and material properties during high-pressure rock failure experiments. It is important to understand the physical quantities that acoustic emission sensors measure, as well as the response of these sensors as a function of frequency. This study calibrates the newly built AE system in the MIT Rock Physics Laboratory using a ball-bouncing system. Full waveforms of multi-bounce events due to ball drops are used to infer the transfer function of lead zirconate titanate (PZT) sensors in high pressure environments. Uncertainty in the sensor transfer functions is quantified using a waveform-based Bayesian approach. The quantification of \textit{in situ} sensor transfer functions makes it possible to apply full waveform analysis for acoustic emissions at high pressures.

BibTeX Citation

@article{ gu-waveform-2020,
  title = { Bayesian waveform-based calibration of high-pressure acoustic emission systems with ball drop measurements },
  author = { C. Gu and U. Mok and Y. M. Marzouk and G. Prieto Gomez and F. Sheibani and J. B. Evans and B. Hager },
  journal = { Geophysical Journal International },
  volume = { 221 },
  pages = { 20-36 },
  month = { 0 },
  year = { 2020 },
  editor = {  },
}

Available Links

• https://arxiv.org/abs/1809.05567

Abstract

We propose a multifidelity dimension reduction method to identify a low-dimensional structure present in many engineering models. The structure of interest arises when functions vary primarily on a low-dimensional subspace of the high-dimensional input space, while varying little along the complementary directions. Our approach builds on the gradient-based methodology of active subspaces, and exploits models of different fidelities to reduce the cost of performing dimension reduction through the computation of the active subspace matrix. We provide a non-asymptotic analysis of the number of gradient evaluations sufficient to achieve a prescribed error in the active subspace matrix, both in expectation and with high probability. We show that the sample complexity depends on a notion of intrinsic dimension of the problem, which can be much smaller than the dimension of the input space. We illustrate the benefits of such a multifidelity dimension reduction approach using numerical experiments with input spaces of up to two thousand dimensions.

BibTeX Citation

@article{ lam-dimred-2020,
  keywords = { Dimension reduction, multifidelity, gradient-based, active subspace, intrinsic dimension, effective rank, matrix Bernstein inequality, control variate },
  title = { Multifidelity dimension reduction via active subspaces },
  author = { R. Lam and O. Zahm and Y. M. Marzouk and K. Willcox },
  journal = { SIAM Journal on Scientific Computing },
  volume = { 42 },
  number = { 2 },
  pages = { A929-A956 },
  month = { 0 },
  year = { 2020 },
  editor = {  },
}

Available Links

• https://doi.org/10.1088/1361-6587/ab4e69

Abstract

We present an analysis which suggests that model selection is a critical ingredient for successful reconstruction of impurity transport coefficient profiles, $D$ and $V$, from experimental data. Determining these quantities is a challenging nonlinear inverse problem. We use synthetic data to show that this problem is ill-posed, and hence $D$ and $V$ are not recommended for use in validation metrics unless the data analysis procedure goes to great lengths to account for the possibility that there are multiple possible solutions. In particular, inferred profiles which are very different from the true ones yield seemingly reasonable goodness-of-fit for synthetic X-ray spectrometer data. We present a Bayesian approach for inferring $D$ and $V$ which provides a rigorous means of selecting the level of complexity of the inferred profiles, thereby enabling successful reconstruction of the profiles.

BibTeX Citation

@article{ art_75,
  title = { On the importance of model selection when inferring impurity transport coefficient profiles },
  editor = {  },
  author = { M. Chilenski and M. Greenwald and Y. M. Marzouk and J. Rice and A. E. White },
  journal = { Plasma Physics and Controlled Fusion },
  doi = { 10.1088/1361-6587/ab4e69 },
  year = { 2019 },
  volume = { 61 },
  month = { 11 },
  pages = { 125012 },
}

Available Links

• https://arxiv.org/abs/1705.04678

Abstract

The development of chemical reaction models aids understanding and prediction in areas ranging from biology to electrochemistry and combustion. A systematic approach to building reaction network models uses observational data not only to estimate unknown parameters, but also to learn model structure. Bayesian inference provides a natural approach to this data-driven construction of models. Yet traditional Bayesian model inference methodologies that numerically evaluate the evidence for each model are often infeasible for nonlinear reaction network inference, as the number of plausible models can be combinatorially large. Alternative approaches based on model-space sampling can enable large-scale network inference, but their realization presents many challenges. In this paper, we present new computational methods that make large-scale nonlinear network inference tractable. First, we exploit the topology of networks describing potential interactions among chemical species to design improved "between-model" proposals for reversible-jump Markov chain Monte Carlo. Second, we introduce a sensitivity-based determination of move types which, when combined with network-aware proposals, yields significant additional gains in sampling performance. These algorithms are demonstrated on inference problems drawn from systems biology, with nonlinear differential equation models of species interactions.

BibTeX Citation

@article{ art_72,
  title = { Exploiting network topology for large-scale inference of nonlinear reaction models },
  keywords = { reaction network, network inference, model selection, Bayesian inference, reversible-jump MCMC. },
  editor = {  },
  author = { N. Galagali and Y. M. Marzouk },
  journal = { Journal of the Royal Society: Interface },
  doi = { https://doi.org/10.1098/rsif.2018.0766 },
  year = { 2019 },
  arxiv = { 1705.04678 },
  volume = { 16 },
  number = { 152 },
  month = { 3 },
  pages = { 20180766 },
}

Available Links

• http://arxiv.org/abs/1510.09088

Abstract

We develop new approximation algorithms and data structures for representing and computing with multivariate functions using the functional tensor-train (FT), a continuous extension of the tensor-train (TT) decomposition. The FT represents functions using a tensor-train ansatz by replacing the three-dimensional TT cores with univariate matrix-valued functions. The main contribution of this paper is a framework to compute the FT that employs adaptive approximations of univariate fibers, and that is not tied to any tensorized discretization. The algorithm can be coupled with any univariate linear or nonlinear approximation procedure. We demonstrate that this approach can generate multivariate function approximations that are several orders of magnitude more accurate, for the same cost, than those based on the conventional approach of compressing the coefficient tensor of a tensor-product basis. Our approach is in the spirit of other continuous computation packages such as Chebfun, and yields an algorithm which requires the computation of “continuous” matrix factorizations such as the LU and QR decompositions of vector-valued functions. To support these developments, we describe continuous versions of an approximate maximum-volume cross approximation algorithm and of a rounding algorithm that re-approximates an FT by one of lower ranks. We demonstrate that our technique improves accuracy and robustness, compared to TT and quantics-TT approaches with fixed parameterizations, of high-dimensional integration, differentiation, and approximation of functions with local features such as discontinuities and other nonlinearities.

BibTeX Citation

@article{ art_70,
  title = { A continuous analogue of the tensor-train decomposition },
  keywords = { tensor decompositions, tensor-train, Chebfun, high-dimensional approximation },
  editor = {  },
  author = { A. Gorodetsky and S. Karaman and Y. M. Marzouk },
  journal = { Computer Methods in Applied Mechanics and Engineering },
  doi = { 10.1016/j.cma.2018.12.015 },
  year = { 2019 },
  arxiv = { 1510.09088 },
  volume = { 347 },
  month = { 0 },
  pages = { 59--84 },
}

Available Links

• http://arxiv.org/abs/1806.03085
• http://papers.neurips.cc/paper/8130-a-stein-variational-newton-method

Abstract

Stein variational gradient descent (SVGD) was recently proposed as a general purpose nonparametric variational inference algorithm [Liu & Wang, NIPS 2016]: it minimizes the Kullback-Leibler divergence between the target distribution and its approximation by implementing a form of functional gradient descent on a reproducing kernel Hilbert space. In this paper, we accelerate and generalize the SVGD algorithm by including second-order information, thereby approximating a Newton-like iteration in function space. We also show how second-order information can lead to more effective choices of kernel. We observe significant computational gains over the original SVGD algorithm in multiple test cases.

BibTeX Citation

@article{ detomasso-svn-2020,
  keywords = { variational inference, Stein's method, transport, nonparametric approximation, gradient descent, Newton method },
  title = { A Stein variational Newton method },
  author = { G. Detomasso and T. Cui and Y. M. Marzouk and R. Scheichl and A. Spantini },
  journal = { Advances in Neural Information Processing Systems (NeurIPS) },
  volume = { 31 },
  month = { 0 },
  year = { 2018 },
  editor = {  },
}

Available Links

• https://arxiv.org/abs/1705.07646

Abstract

We study Bayesian inference methods for solving linear inverse problems, focusing on hierarchical formulations where the prior or the likelihood function depend on unspecified hyperparameters. In practice, these hyperparameters are often determined via an empirical Bayesian method that maximizes the marginal likelihood function, i.e., the probability density of the data conditional on the hyperparameters. Evaluating the marginal likelihood, however, is computationally challenging for large-scale problems. In this work, we present a method to approximately evaluate marginal likelihood functions, based on a low-rank approximation of the update from the prior covariance to the posterior covariance. We show that this approximation is optimal in a minimax sense. Moreover, we provide an efficient algorithm to implement the proposed method, based on a combination of the randomized SVD and a spectral approximation method to compute square roots of the prior covariance matrix. Several numerical examples demonstrate good performance of the proposed method.

BibTeX Citation

@article{ art_65,
  title = { An approximate empirical Bayesian method for large-scale linear-Gaussian inverse problems },
  editor = {  },
  author = { Q. Zhou and W. Liu and J. Li and Y. M. Marzouk },
  journal = { Inverse Problems },
  doi = { 10.1088/1361-6420/aac287 },
  year = { 2018 },
  arxiv = { 1705.07646 },
  volume = { 34 },
  number = { 9 },
  month = { 0 },
  pages = { 095001 },
}

Available Links

• https://doi.org/10.1016/j.patcog.2018.03.022

Abstract

Novelty detection is an important task in a variety of applications such as object recognition, defect localization, medical diagnostics, and event detection. The objective of novelty detection is to distinguish one class, for which data are available, from all other possible classes when there is insufficient information to build an explicit model for the latter. The data from the observed class are usually represented in terms of certain features which can be modeled as random variables (RV). An important challenge for novelty detection in multivariate problems is characterizing the statistical dependencies among these RVs. Failure to consider these dependencies may lead to inaccurate predictions, usually in the form of high false positive rates. In this study, we propose conditional classifiers as a new approach for novelty detection that is capable of accounting for statistical dependencies of the relevant RVs without simplifying assumptions. To implement the proposed idea, we use Gaussian mixture models (GMM) along with forward stage-wise additive modeling and boosting methods to learn the conditional densities of RVs that represent our observed data. The resulting model, which is called a boosted conditional GMM, is then used as a basis for classification. To test the performance of the proposed method, we apply it to a realistic application problem for analyzing sensor networks and compare the results to those obtained with a one-class support vector machine.

BibTeX Citation

@article{ art_63,
  title = { Conditional classifiers and boosted conditional Gaussian mixture models for novelty detection },
  keywords = { Novelty detection, mixture models, graphical models, conditional dependence, conditional density, additive modeling, boosting, false positive },
  editor = {  },
  author = { R. Mohammadi-Ghazi and Y. M. Marzouk and O. Büyüköztürk },
  journal = { Pattern Recognition },
  doi = { 10.1016/j.patcog.2018.03.022 },
  year = { 2018 },
  volume = { 81 },
  month = { 0 },
  pages = { 601--614 },
}

Available Links

• https://doi.org/10.1063/1.4997251

Abstract

Recent attempts to measure impurity transport in Alcator C-Mod using an x-ray imaging crystal spectrometer and laser blow-off impurity injector have failed to yield unique reconstructions of the transport coefficient profiles. This paper presents a fast, linearized model which was constructed to estimate diagnostic requirements for impurity transport experiments. The analysis shows that the spectroscopic diagnostics on Alcator C-Mod should be capable of inferring simple profiles of impurity diffusion $D_Z$ and convection $V_Z$ accurate to better than ±10% uncertainty, suggesting that the failure to infer unique $D_Z$, $V_Z$ from experimental data is attributable to an inadequate analysis procedure rather than the result of insufficient diagnostics. Furthermore, the analysis reveals that spatial resolution is more important than temporal resolution for typical diagnostic sampling rates and noise levels. This approach can be adapted to design and verify diagnostics for transport experiments on any magnetic confinement device.

BibTeX Citation

@article{ art_61,
  title = { Efficient design and verification of diagnostics for impurity transport experiments },
  editor = {  },
  author = { M. Chilenski and M. Greenwald and Y. M. Marzouk and J. Rice and A. E. White },
  journal = { Review of Scientific Instruments },
  doi = { 10.1063/1.4997251 },
  year = { 2018 },
  volume = { 89 },
  month = { 1 },
  pages = { 013504 },
}

Available Links

• https://doi.org/10.1177/0278364917753994

Abstract

Motion planning and control problems are embedded and essential in almost all robotics applications. These problems are often formulated as stochastic optimal control problems and solved using dynamic programming algorithms. Unfortunately, most existing algorithms that guarantee convergence to optimal solutions suffer from the curse of dimensionality: the run time of the algorithm grows exponentially with the dimension of the state space of the system. We propose novel dynamic programming algorithms that alleviate the curse of dimensionality in problems that exhibit certain low-rank structure. The proposed algorithms are based on continuous tensor decompositions recently developed by the authors. Essentially, the algorithms represent high-dimensional functions (e.g., the value function) in a compressed format, and directly perform dynamic programming computations (e.g., value iteration, policy iteration) in this format. Under certain technical assumptions, the new algorithms guarantee convergence towards optimal solutions with arbitrary precision. Furthermore, the run times of the new algorithms scale polynomially with the state dimension and polynomially with the ranks of the value function. This approach realizes substantial computational savings in “compressible” problem instances, where value functions admit low-rank approximations. We demonstrate the new algorithms in a wide range of problems, including a simulated six-dimensional agile quadcopter maneuvering example and a seven-dimensional aircraft perching example. In some of these examples, we estimate computational savings of up to 10 orders of magnitude over standard value iteration algorithms. We further demonstrate the algorithms running in real time on board a quadcopter during a flight experiment under motion capture.

BibTeX Citation

@article{ art_59,
  title = { High-dimensional stochastic optimal control using continuous tensor decompositions },
  keywords = { stochastic optimal control, motion planning, dynamic programming, tensor decompositions },
  editor = {  },
  author = { A. Gorodetsky and S. Karaman and Y. M. Marzouk },
  journal = { International Journal of Robotics Research },
  doi = { 10.1177/0278364917753994 },
  year = { 2018 },
  arxiv = { arXiv:1611.04706 },
  volume = { 37 },
  number = { 2-3 },
  month = { 3 },
  pages = { 340--377 },
}

Available Links

• http://jmlr.org/papers/v19/17-747.html

Abstract

Integration against an intractable probability measure is among the fundamental challenges of statistical inference, particularly in the Bayesian setting. A principled approach to this problem seeks a deterministic coupling of the measure of interest with a tractable ``reference'' measure (e.g., a standard Gaussian). This coupling is induced by a transport map, and enables direct simulation from the desired measure simply by evaluating the transport map at samples from the reference. Yet characterizing such a map---e.g., representing and evaluating it---grows challenging in high dimensions. The central contribution of this paper is to establish a link between the Markov properties of the target measure and the existence of certain low-dimensional couplings, induced by transport maps that are sparse or decomposable. Our analysis not only facilitates the construction of couplings in high-dimensional settings, but also suggests new inference methodologies. For instance, in the context of nonlinear and non-Gaussian state space models, we describe new variational algorithms for online filtering, smoothing, and parameter estimation. These algorithms implicitly characterize---via a transport map---the full posterior distribution of the sequential inference problem using local operations only incrementally more complex than regular filtering, while avoiding importance sampling or resampling.

BibTeX Citation

@article{ art_57,
  keywords = { transport map, rearrangement, Bayesian inference, variational inference, graphical model, Markov random field, sparsity, Kalman recursions, filtering, smoothing, joint parameter-state estimation, state-space model },
  title = { Inference via low-dimensional couplings },
  author = { A. Spantini and D. Bigoni and Y. M. Marzouk },
  journal = { The Journal of Machine Learning Research },
  volume = { 19 },
  number = { 66 },
  pages = { 1--71 },
  month = { 0 },
  year = { 2018 },
  editor = {  },
}

Available Links

• https://arc.aiaa.org/doi/abs/10.2514/1.J055877

Abstract

Many simulation tools for airfoil analysis and design are based on an integral approximation of the boundary layer. This approximate formulation cannot resolve the full dynamics of boundary-layer flows and hence requires additional models to account for unresolved effects. This paper introduces a new, data-driven, probabilistic model of these unresolved effects for the incompressible and laminar regime. To construct this model, methods from supervised learning have been applied to a dataset containing over 1550 airfoils. The result is a model that 1) is based on a large dataset of realistic airfoil configurations, and 2) quantifies the model inadequacy associated with the use of an approximate boundary-layer formulation. A stochastic version of the airfoil design tool XFOIL has also been created by replacing its original boundary-layer model with the probabilistic model developed here. This stochastic version of XFOIL has been applied to compute the drag polars of two airfoils at low Reynolds numbers and the results are compared with experimental data.

BibTeX Citation

@article{ art_53,
  title = { Data-driven integral boundary layer modeling for airfoil performance prediction in the laminar regime },
  editor = {  },
  author = { A. Marques and Q. Wang and Y. M. Marzouk },
  journal = { AIAA Journal },
  doi = { 10.2514/1.J055877 },
  year = { 2018 },
  volume = { 56 },
  number = { 2 },
  month = { 0 },
  pages = { 482--496 },
}

Available Links

• http://papers.nips.cc/paper/6830-beyond-normality-learning-sparse-probabilistic-graphical-models-in-the-non-gaussian-setting

Abstract

We present an algorithm to identify sparse dependence structure in continuous and non-Gaussian probability distributions, given a corresponding set of data. The conditional independence structure of an arbitrary distribution can be represented as an undirected graph (or Markov random field), but most algorithms for learning this structure are restricted to the discrete or Gaussian cases. Our new approach allows for more realistic and accurate descriptions of the distribution in question, and in turn better estimates of its sparse Markov structure. Sparsity in the graph is of interest as it can accelerate inference, improve sampling methods, and reveal important dependencies between variables. The algorithm relies on exploiting the connection between the sparsity of the graph and the sparsity of transport maps, which deterministically couple one probability measure to another.

BibTeX Citation

@article{ art_51,
  title = { Beyond normality: Learning sparse probabilistic graphical models in the non-Gaussian setting },
  author = { R. Morrison and R. Baptista and Y. M. Marzouk },
  journal = { Advances in Neural Information Processing Systems (NIPS) },
  volume = { 30 },
  month = { 0 },
  year = { 2017 },
  editor = {  },
}

Available Links

• http://arxiv.org/abs/1607.01881

Abstract

We propose optimal dimensionality reduction techniques for the solution of goal-oriented linear-Gaussian inverse problems, where the quantity of interest (QoI) is a function of the inversion parameters. These approximations are suitable for large-scale applications. In particular, we study the approximation of the posterior covariance of the QoI as a low-rank negative update of its prior covariance, and prove optimality of this update with respect to the natural geodesic distance on the manifold of symmetric positive definite matrices. Assuming exact knowledge of the posterior mean of the QoI, the optimality results extend to optimality in distribution with respect to the Kullback-Leibler divergence and the Hellinger distance between the associated distributions. We also propose approximation of the posterior mean of the QoI as a low-rank linear function of the data, and prove optimality of this approximation with respect to a weighted Bayes risk. Both of these optimal approximations avoid the explicit computation of the full posterior distribution of the parameters and instead focus on directions that are well informed by the data and relevant to the QoI. These directions stem from a balance among all the components of the goal-oriented inverse problem: prior information, forward model, measurement noise, and ultimate goals. We illustrate the theory using a high-dimensional inverse problem in heat transfer.

BibTeX Citation

@article{ art_49,
  title = { Goal-oriented optimal approximations of Bayesian linear inverse problems },
  keywords = { inverse problems, goal--oriented, Bayesian inference, low-rank approximation, covariance approximation, Riemannian metric, geodesic distance, posterior mean approximation, Bayes risk, optimality },
  editor = {  },
  author = { A. Spantini and T. Cui and K. Willcox and L. Tenorio and Y. M. Marzouk },
  journal = { SIAM Journal on Scientific Computing },
  doi = { 10.1137/16M1082123 },
  year = { 2017 },
  arxiv = { 1607.01881 },
  volume = { 39 },
  number = { 5 },
  month = { 0 },
  pages = { S167--S196 },
}

Available Links

• http://dx.doi.org/10.1137/15M1032478

Abstract

In many inverse problems, model parameters cannot be precisely determined from observational data. Bayesian inference provides a mechanism for capturing the resulting parameter uncertainty, but typically at a high computational cost. This work introduces a multiscale decomposition that exploits conditional independence across scales, when present in certain classes of inverse problems, to decouple Bayesian inference into two stages: (1) a computationally tractable coarse-scale inference problem; and (2) a mapping of the low-dimensional coarse-scale posterior distribution into the original high-dimensional parameter space. This decomposition relies on a characterization of the non-Gaussian joint distribution of coarse- and fine-scale quantities via optimal transport maps. We demonstrate our approach on a sequence of inverse problems arising in subsurface flow, using the multiscale finite element method to discretize the steady state pressure equation. We compare the multiscale strategy with full-dimensional Markov chain Monte Carlo on a problem of moderate dimension (100 parameters) and then use it to infer a conductivity field described by over 10,000 parameters.

BibTeX Citation

@article{ parno-transport-juq-2016,
  keywords = { Bayesian inference; inverse problems; multiscale modeling; multiscale finite element method; optimal transportation; Markov chain Monte Carlo },
  doi = { 10.1137/15M1032478 },
  title = { A multiscale strategy for Bayesian inference using transport maps },
  author = { M. Parno and T. Moselhy and Y. M. Marzouk },
  journal = { SIAM/ASA Journal on Uncertainty Quantification },
  volume = { 4 },
  number = { 1 },
  pages = { 1160--1190 },
  month = { 10 },
  year = { 2016 },
  editor = {  },
}

Available Links

• http://arxiv.org/abs/1510.06053

Abstract

Two major bottlenecks to the solution of large-scale Bayesian inverse problems are the scaling of posterior sampling algorithms to high dimensional parameter spaces and the computational cost of forward model evaluations. Yet incomplete or noisy data, the state variation and parameter dependence of the forward model, and correlations in the prior collectively provide useful structure that can be exploited for dimension reduction in this setting---both in the parameter space of the inverse problem and in the state space of the forward model. To this end, we show how to jointly construct low-dimensional subspaces of the parameter space and the state space in order to accelerate the Bayesian solution of the inverse problem. As a byproduct of state dimension reduction, we also show how to identify low-dimensional subspaces of the data in problems with high-dimensional observations. These subspaces enable approximation of the posterior as a product of two factors: (i) a projection of the posterior onto a low-dimensional parameter subspace, wherein the original likelihood is replaced by an approximation involving a reduced model; and (ii) the marginal prior distribution on the high-dimensional complement of the parameter subspace. We present and compare several strategies for constructing these subspaces using only a limited number of forward and adjoint model simulations. The resulting posterior approximations can rapidly be characterized using standard sampling techniques, e.g., Markov chain Monte Carlo. Two numerical examples demonstrate the accuracy and efficiency of our approach: inversion of an integral equation in atmospheric remote sensing, where the data dimension is very high; and the inference of a heterogeneous transmissivity field in a groundwater system, which involves a partial differential equation forward model with high dimensional state and parameters.

BibTeX Citation

@article{ art_44,
  title = { Scalable posterior approximations for large-scale Bayesian inverse problems via likelihood-informed parameter and state reduction },
  keywords = { Inverse problems, Bayesian inference, dimension reduction, model reduction, low-rank approximation, Markov chain Monte Carlo },
  editor = {  },
  author = { T. Cui and Y. M. Marzouk and K. Willcox },
  journal = { Journal of Computational Physics },
  doi = { doi:10.1016/j.jcp.2016.03.055 },
  year = { 2016 },
  arxiv = { 1510.06053 },
  volume = { 315 },
  month = { 6 },
  pages = { 363--387 },
}

Available Links

• http://arxiv.org/abs/1508.06452

Abstract

A priori dimension reduction is a widely adopted technique for reducing the computational complexity of stationary inverse problems. In this setting, the solution of an inverse problem is parameterized by a low-dimensional basis that is often obtained from the truncated Karhunen-Loève expansion of the prior distribution. For high-dimensional inverse problems equipped with smoothing priors, this technique can lead to drastic reductions in parameter dimension and significant computational savings. In this paper, we extend the concept of a priori dimension reduction to non-stationary inverse problems, in which the goal is to sequentially infer the state of a dynamical system. Our approach proceeds in an offline-online fashion. We first identify a low-dimensional subspace in the state space before solving the inverse problem (the offline phase), using either the method of "snapshots" or regularized covariance estimation. Then this subspace is used to reduce the computational complexity of various filtering algorithms—including the Kalman filter, extended Kalman filter, and ensemble Kalman filter—within a novel subspace-constrained Bayesian prediction-and-update procedure (the online phase). We demonstrate the performance of our new dimension reduction approach on various numerical examples. In some test cases, our approach reduces the dimensionality of the original problem by orders of magnitude and yields up to two orders of magnitude in computational savings.

BibTeX Citation

@article{ art_42,
  title = { On dimension reduction in Gaussian filters },
  editor = {  },
  author = { A. Solonen and T. Cui and J. Hakkarainen and Y. M. Marzouk },
  journal = { Inverse Problems },
  doi = { 10.1088/0266-5611/32/4/045003 },
  year = { 2016 },
  arxiv = { 1508.06452 },
  volume = { 32 },
  number = { 4 },
  month = { 3 },
  pages = { 045003 },
}

Available Links

• http://arxiv.org/abs/1411.3688

Abstract

Many Bayesian inference problems require exploring the posterior distribution of high-dimensional parameters that represent the discretization of an underlying function. This work introduces a family of Markov chain Monte Carlo (MCMC) samplers that can adapt to the particular structure of a posterior distribution over functions. Two distinct lines of research intersect in the methods developed here. First, we introduce a general class of operator-weighted proposal distributions that are well defined on function space, such that the performance of the resulting MCMC samplers is independent of the discretization of the function. Second, by exploiting local Hessian information and any associated low-dimensional structure in the change from prior to posterior distributions, we develop an inhomogeneous discretization scheme for the Langevin stochastic differential equation that yields operator-weighted proposals adapted to the non-Gaussian structure of the posterior. The resulting dimension-independent, likelihood-informed (DILI) MCMC samplers may be useful for a large class of high-dimensional problems where the target probability measure has a density with respect to a Gaussian reference measure. Two nonlinear inverse problems are used to demonstrate the efficiency of these DILI samplers: an elliptic PDE coefficient inverse problem and path reconstruction in a conditioned diffusion.

BibTeX Citation

@article{ art_38,
  title = { Dimension-independent likelihood-informed MCMC },
  keywords = { Markov chain Monte Carlo, likelihood-informed subspace, infinite-dimensional inverse problems, Langevin SDE, conditioned diffusion },
  editor = {  },
  author = { T. Cui and K. J. H. Law and Y. M. Marzouk },
  journal = { Journal of Computational Physics },
  doi = { doi:10.1016/j.jcp.2015.10.008 },
  year = { 2016 },
  arxiv = { 1411.3688 },
  volume = { 304 },
  month = { 1 },
  pages = { 109--137 },
}

Available Links

• http://arxiv.org/abs/1407.3463

Abstract

In the Bayesian approach to inverse problems, data are often informative, relative to the prior, only on a low-dimensional subspace of the parameter space. Significant computational savings can be achieved by using this subspace to characterize and approximate the posterior distribution of the parameters. We first investigate approximation of the posterior covariance matrix as a low-rank update of the prior covariance matrix. We prove optimality of a particular update, based on the leading eigendirections of the matrix pencil defined by the Hessian of the negative log-likelihood and the prior precision, for a broad class of loss functions. This class includes the F\"{o}rstner metric for symmetric positive definite matrices, as well as the Kullback-Leibler divergence and the Hellinger distance between the associated distributions. We also propose two fast approximations of the posterior mean and prove their optimality with respect to a weighted Bayes risk under squared-error loss. These approximations are deployed in an offline-online manner, where a more costly but data-independent offline calculation is followed by fast online evaluations. As a result, these approximations are particularly useful when repeated posterior mean evaluations are required for multiple data sets. We demonstrate our theoretical results with several numerical examples, including high-dimensional X-ray tomography and an inverse heat conduction problem. In both of these examples, the intrinsic low-dimensional structure of the inference problem can be exploited while producing results that are essentially indistinguishable from solutions computed in the full space.

BibTeX Citation

@article{ art_35,
  title = { Optimal low-rank approximations of Bayesian linear inverse problems },
  keywords = { inverse problems, Bayesian inference, low-rank approximation, covariance approximation, Förstner-Moonen metric, posterior mean approximation, Bayes risk, optimality },
  editor = {  },
  author = { A. Spantini and A. Solonen and T. Cui and J. Martin and L. Tenorio and Y. M. Marzouk },
  journal = { SIAM Journal on Scientific Computing },
  doi = { 10.1137/140977308 },
  year = { 2015 },
  arxiv = { 1407.3463 },
  volume = { 37 },
  number = { 6 },
  month = { 0 },
  pages = { A2451--A2487 },
}

Available Links

• http://dx.doi.org/10.1002/nme.4928

Abstract

We introduce a novel numerical approach to Bayesian parameter estimation in partial differential equations. The main idea is to translate the equation into a state-discrete dynamic Bayesian network with the discretization of cellular probabilistic automata. There exists a vast pool of inference algorithms in the probabilistic graphical models framework which can be applied to the network. In particular, we reformulate the parameter estimation problem as a filtering problem, discuss requirements for our specific setup, and apply the Boyen-Koller algorithm. To demonstrate our ideas, the scheme is applied to the problem of arsenate advection and adsorption in a water pipe: from measurements of the concentration of dissolved arsenate at the outflow boundary condition, we infer the strength of an arsenate source at the inflow boundary condition.

BibTeX Citation

@article{ art_33,
  title = { A new network approach to Bayesian inference in partial differential equations },
  keywords = { partial differential equations, dynamic Bayesian networks, Boyen--Koller algorithm, cellular probabilistic automata },
  editor = {  },
  author = { D. Kohler and Y. M. Marzouk and J. Müller and U. Wever },
  journal = { International Journal for Numerical Methods in Engineering },
  doi = { 10.1002/nme.4928 },
  year = { 2015 },
  volume = { 104 },
  number = { 5 },
  month = { 11 },
  pages = { 313--329 },
}

Available Links

• http://stacks.iop.org/0965-0393/23/015009

Abstract

We demonstrate that, by observing the behavior of a film deposited on a substrate, certain features of the substrate may be inferred with quantified uncertainty using Bayesian methods. We carry out this demonstration on an illustrative film/substrate model, where the substrate is a Gaussian random field and the film is a two-component mixture that obeys the Cahn-Hilliard equation. We construct a stochastic reduced order model to describe the film/substrate interaction and use it to infer substrate properties from film behavior. This quantitative inference strategy may be adapted to other film/substrate systems.

BibTeX Citation

@article{ art_31,
  title = { Bayesian inference of substrate properties from film behavior },
  editor = {  },
  author = { R. Aggarwal and M. Demkowicz and Y. M. Marzouk },
  journal = { Modelling and Simulation in Materials Science and Engineering },
  doi = { http://dx.doi.org/10.1088/0965-0393/23/1/015009 },
  year = { 2015 },
  volume = { 23 },
  month = { 0 },
  pages = { 015009 },
}

Available Links

• http://arxiv.org/abs/1403.4290

Abstract

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.

BibTeX Citation

@article{ art_29,
  title = { Data-driven model reduction for the Bayesian solution of inverse problems },
  keywords = { model reduction; inverse problem; adaptive Markov chain Monte Carlo; approximate Bayesian inference },
  editor = {  },
  author = { T. Cui and Y. M. Marzouk and K. Willcox },
  journal = { International Journal for Numerical Methods in Engineering },
  doi = { doi:10.1002/nme.4748 },
  year = { 2015 },
  arxiv = { 1403.4290 },
  volume = { 102 },
  number = { 5 },
  month = { 0 },
  pages = { 966--990 },
}

Available Links

• http://arxiv.org/abs/1402.2845

Abstract

Surrogate models for computational simulations are input-output approximations that allow computationally intensive analyses, such as uncertainty propagation and inference, to be performed efficiently. When a simulation output does not depend smoothly on its inputs, the error and convergence rate of many approximation methods deteriorate substantially. This paper details a method for efficiently localizing discontinuities in the input parameter domain, so that the model output can be approximated as a piecewise smooth function. The approach comprises an initialization phase, which uses polynomial annihilation to assign function values to different regions and thus seed an automated labeling procedure, followed by a refinement phase that adaptively updates a kernel support vector machine representation of the separating surface via active learning. The overall approach avoids structured grids and exploits any available simplicity in the geometry of the separating surface, thus reducing the number of model evaluations required to localize the discontinuity. The method is illustrated on examples of up to eleven dimensions, including algebraic models and ODE/PDE systems, and demonstrates improved scaling and efficiency over other discontinuity localization approaches.

BibTeX Citation

@article{ art_27,
  title = { Efficient localization of discontinuities in complex computational simulations },
  keywords = { discontinuity detection, polynomial annihilation, function approximation, support vector machines, active learning, uncertainty quantification },
  editor = {  },
  author = { A. Gorodetsky and Y. M. Marzouk },
  journal = { SIAM Journal on Scientific Computing },
  year = { 2014 },
  arxiv = { 1402.2845 },
  volume = { 36 },
  number = { 6 },
  month = { 0 },
  pages = { A2584--A2610 },
}

Available Links

• http://arxiv.org/abs/1212.2228

Abstract

Optimal experimental design (OED) seeks experiments expected to yield the most useful data for some purpose. In practical circumstances where experiments are time-consuming or resource-intensive, OED can yield enormous savings. We pursue OED for nonlinear systems from a Bayesian perspective, with the goal of choosing experiments that are optimal for parameter inference. Our objective in this context is the expected information gain in model parameters, which in general can only be estimated using Monte Carlo methods. Maximizing this objective thus becomes a stochastic optimization problem. This paper develops gradient-based stochastic optimization methods for the design of experiments on a continuous parameter space. Given a Monte Carlo estimator of expected information gain, we use infinitesimal perturbation analysis to derive gradients of this estimator. We are then able to formulate two gradient-based stochastic optimization approaches: (i) Robbins-Monro stochastic approximation, and (ii) sample average approximation combined with a deterministic quasi-Newton method. A polynomial chaos approximation of the forward model accelerates objective and gradient evaluations in both cases. We discuss the implementation of these optimization methods, then conduct an empirical comparison of their performance. To demonstrate design in a nonlinear setting with partial differential equation forward models, we use the problem of sensor placement for source inversion. Numerical results yield useful guidelines on the choice of algorithm and sample sizes, assess the impact of estimator bias, and quantify tradeoffs of computational cost versus solution quality and robustness.

BibTeX Citation

@article{ art_25,
  title = { Gradient-based stochastic optimization methods in Bayesian experimental design },
  editor = {  },
  author = { X. Huan and Y. M. Marzouk },
  journal = { International Journal for Uncertainty Quantification },
  doi = { 10.1615/Int.J.UncertaintyQuantification.2014006730 },
  year = { 2014 },
  arxiv = { http://arxiv.org/abs/1212.2228 },
  volume = { 4 },
  number = { 6 },
  month = { 0 },
  pages = { 479--510 },
  comments = { International Journal for Uncertainty Quantification, in press },
}

Available Links

• http://arxiv.org/abs/1209.1406

Abstract

Polynomial approximations of computationally intensive models are central to uncertainty quantification. This paper describes an adaptive method for non-intrusive pseudospectral approximation, based on Smolyak's algorithm with generalized sparse grids. We rigorously analyze and extend the non-adaptive method proposed in [Constantine et al. 2012], and compare it to a common alternative approach for using sparse grids to construct polynomial approximations, direct quadrature. Analysis of direct quadrature shows that O(1) errors are an intrinsic property of some configurations of the method, as a consequence of internal aliasing. We provide precise conditions, based on the chosen polynomial basis and quadrature rules, under which this aliasing error occurs. We then establish theoretical results on the accuracy of Smolyak pseudospectral approximation, and show that the Smolyak approximation avoids internal aliasing and makes far more effective use of sparse function evaluations. These results are applicable to broad choices of quadrature rule and generalized sparse grids. Exploiting this flexibility, we introduce a greedy heuristic for adaptive refinement of the pseudospectral approximation. We numerically demonstrate convergence of the algorithm on the Genz test functions, and illustrate the accuracy and efficiency of the adaptive approach on a realistic chemical kinetics problem.

BibTeX Citation

@article{ art_23,
  title = { Adaptive Smolyak pseudospectral approximations },
  keywords = { Smolyak algorithms, sparse grids, orthogonal polynomials, pseudospectral approximation, approximation theory, uncertainty quantification },
  editor = {  },
  author = { P. Conrad and Y. M. Marzouk },
  journal = { SIAM Journal on Scientific Computing },
  doi = { 10.1137/120890715 },
  year = { 2013 },
  arxiv = { http://arxiv.org/abs/1209.1406 },
  volume = { 35 },
  number = { 6 },
  month = { 0 },
  pages = { A2643--A2670 },
}

Available Links

• http://arxiv.org/abs/1108.4146

Abstract

The optimal selection of experimental conditions is essential to maximizing the value of data for inference and prediction, particularly in situations where experiments are time-consuming and expensive to conduct. We propose a general mathematical framework and an algorithmic approach for optimal experimental design with nonlinear simulation-based models; in particular, we focus on finding sets of experiments that provide the most information about targeted sets of parameters. Our framework employs a Bayesian statistical setting, which provides a foundation for inference from noisy, indirect, and incomplete data, and a natural mechanism for incorporating heterogeneous sources of information. An objective function is constructed from information theoretic measures, reflecting expected information gain from proposed combinations of experiments. Polynomial chaos approximations and a two-stage Monte Carlo sampling method are used to evaluate the expected information gain. Stochastic approximation algorithms are then used to make optimization feasible in computationally intensive and high-dimensional settings. These algorithms are demonstrated on model problems and on nonlinear parameter inference problems arising in detailed combustion kinetics.

BibTeX Citation

@article{ art_21,
  title = { Simulation-based optimal Bayesian experimental design for nonlinear systems },
  keywords = { Uncertainty quantification; Bayesian inference; Optimal experimental design; Nonlinear experimental design; Stochastic approximation; Shannon information; Chemical kinetics },
  editor = {  },
  author = { X. Huan and Y. M. Marzouk },
  journal = { Journal of Computational Physics },
  doi = { 10.1016/j.jcp.2012.08.013 },
  year = { 2013 },
  arxiv = { 1108.4146 },
  volume = { 232 },
  number = { 1 },
  month = { 1 },
  pages = { 288--317 },
  comments = { test },
}

Available Links

• http://dx.doi.org/10.1016/j.advwatres.2012.04.009

Abstract

BibTeX Citation

@article{ art_19,
  title = { Bayesian reconstruction of binary media with unresolved fine-scale spatial structures },
  keywords = { Upscaling; Binary media; Bayesian technique; Multiscale inference },
  editor = {  },
  author = { J. Ray and S. McKenna and B. van Bloemen Waanders and Y. M. Marzouk },
  journal = { Advances in Water Resources },
  doi = { 10.1016/j.advwatres.2012.04.009 },
  year = { 2012 },
  volume = { 44 },
  month = { 8 },
  pages = { 1--19 },
}

Available Links

• http://dx.doi.org/10.1016/j.jcp.2011.10.031

Abstract

A critical problem in accurately estimating uncertainty in model predictions is the lack of details in the literature on the correlation (or full joint distribution) of uncertain model parameters. In this paper we describe a framework and a class of algorithms for analyzing such “missing data” problems in the setting of Bayesian statistics. The analysis focuses on the family of posterior distributions consistent with given statistics (e.g. nominal values, confidence intervals). The combining of consistent distributions is addressed via techniques from the opinion pooling literature. The developed approach allows subsequent propagation of uncertainty in model inputs consistent with reported statistics, in the absence of data.

BibTeX Citation

@article{ art_17,
  title = { Data-free inference of the joint distribution of uncertain model parameters },
  keywords = { Uncertainty quantification; Bayesian statistics; Missing information },
  editor = {  },
  author = { R. D. Berry and H. N. Najm and B. J. Debusschere and Y. M. Marzouk and H. Adalsteinsson },
  journal = { Journal of Computational Physics },
  doi = { 10.1016/j.jcp.2011.10.031 },
  year = { 2012 },
  volume = { 231 },
  number = { 5 },
  month = { 0 },
  pages = { 2180-2198 },
}

Available Links

• http://dx.doi.org/10.1016/j.advwatres.2011.02.011

Abstract

Truncated Gaussian fields provide a flexible model for defining binary media with dispersed (as opposed to layered) inclusions. General properties of excursion sets on these truncated fields are coupled with a distance-based upscaling algorithm and approximations of point process theory to develop an estimation approach for effective conductivity in two-dimensions. Estimation of effective conductivity is derived directly from knowledge of the kernel size used to create the multiGaussian field, defined as the full-width at half maximum (FWHM), the truncation threshold and conductance values of the two modes. Therefore, instantiation of the multiGaussian field is not necessary for estimation of the effective conductance. The critical component of the effective medium approximation developed here is the mean distance between high conductivity inclusions. This mean distance is characterized as a function of the FWHM, the truncation threshold and the ratio of the two modal conductivities. Sensitivity of the resulting effective conductivity to this mean distance is examined for two levels of contrast in the modal conductances and different FWHM sizes. Results demonstrate that the FWHM is a robust measure of mean travel distance in the background medium. The resulting effective conductivities are accurate when compared to numerical results and results obtained from effective media theory, distance-based upscaling and numerical simulation.

BibTeX Citation

@article{ art_15,
  title = { Truncated multi-Gaussian fields and effective conductance of binary media },
  keywords = { Upscaling; Binary media; Effective conductivity },
  editor = {  },
  author = { S. McKenna and J. Ray and Y. M. Marzouk and B. van Bloemen Waanders },
  journal = { Advances in Water Resources },
  doi = { 10.1016/j.advwatres.2011.02.011 },
  year = { 2011 },
  volume = { 34 },
  number = { 5 },
  month = { 5 },
  pages = { 617-626 },
}

Available Links

• http://dx.doi.org/10.1137/10080590X

Abstract

This work focuses on characterizing the integral features of atomic diffusion in Ni/Al nanolaminates based on molecular dynamics (MD) computations. Attention is focused on the simplified problem of extracting the diffusivity, D, in an isothermal system at high temperature. To this end, a mixing measure theory is developed that relies on analyzing the moments of the cumulative distribution functions (CDFs) of the constituents. The mixing measures obtained from replica simulations are exploited in a Bayesian inference framework, based on contrasting these measures with corresponding moments of a dimensionless concentration evolving according to a Fickian process. The noise inherent in the MD simulations is described as a Gaussian process, and this hypothesis is verified both a priori and using a posterior predictive check. Computed values of D for an initially unmixed system rapidly heated to 1500 K are found to be consistent with experimental correlation for diffusion of Ni into molten Al. On the contrary, large discrepancies with experimental predictions are observed when D is estimated based on large-time mean-square displacement (MSD) analysis, and when it is evaluated using the Arrhenius correlation calibrated against experimental measurements of self-propagating front velocities. Implications are finally drawn regarding extension of the present work and potential refinement of continuum modeling approaches.

BibTeX Citation

@article{ art_13,
  title = { Bayesian inference of atomic diffusivity in a binary Ni/Al system based on molecular dynamics },
  editor = {  },
  author = { F. Rizzi and M. Salloum and Y. M. Marzouk and R. Xu and M. L. Falk and T. P. Weihs and G. Fritz and O. M. Knio },
  journal = { SIAM Multiscale Modeling and Simulation },
  doi = { 10.1137/10080590X },
  year = { 2011 },
  volume = { 9 },
  number = { 1 },
  month = { 0 },
  pages = { 486-512 },
}

Available Links

• http://dx.doi.org/10.1137/080726872

Abstract

We study the convergence characteristics of two algebraic kernels used in vortex calculations: the Rosenhead--Moore kernel, which is a low-order kernel, and the Winckelmans--Leonard kernel, which is a high-order kernel. To facilitate the study, a method of evaluating particle-cluster interactions is introduced for the Winckelmans--Leonard kernel. The method is based on Taylor series expansion in Cartesian coordinates, as initially proposed by Lindsay and Krasny [J. Comput. Phys., 172 (2001), pp. 879--907] for the Rosenhead--Moore kernel. A recurrence relation for the Taylor coefficients of the Winckelmans--Leonard kernel is derived by separating the kernel into two parts, and an error estimate is obtained to ensure adaptive error control. The recurrence relation is incorporated into a tree-code to evaluate vorticity-induced velocity. Next, comparison of convergence is made while utilizing the tree-code. Both algebraic kernels lead to convergence, but the Winckelmans--Leonard kernel exhibits a superior convergence rate. The combined desingularization and discretization error from the Winckelmans--Leonard kernel is an order of magnitude smaller than that from the Rosenhead--Moore kernel at a typical resolution. Simulations of vortex rings are performed using the two algebraic kernels in order to compare their performance in a practical setting. In particular, numerical simulations of the side-by-side collision of two identical vortex rings suggest that the three-dimensional evolution of vorticity at finite resolution can be greatly affected by the choice of the kernel. We find that the Winckelmans--Leonard kernel is able to perform the same task with a much smaller number of vortex elements than the Rosenhead--Moore kernel, greatly reducing the overall computational cost.

BibTeX Citation

@article{ art_11,
  title = { Convergence characteristics and computational cost of two algebraic kernels in vortex methods with a tree-code algorithm },
  editor = {  },
  author = { D. Wee and Y. M. Marzouk and F. Schlegel and A. F. Ghoniem },
  journal = { SIAM Journal on Scientific Computing },
  doi = { 10.1137/080726872 },
  year = { 2009 },
  volume = { 31 },
  number = { 4 },
  month = { 0 },
  pages = { 2510-2527 },
}

Available Links

• http://docs.lib.purdue.edu/prism/16/

Abstract

We present an efficient numerical strategy for the Bayesian solution of inverse problems. Stochastic collocation methods, based on generalized polynomial chaos (gPC), are used to construct a polynomial approximation of the forward solution over the support of the prior distribution. This approximation then defines a surrogate posterior probability density that can be evaluated repeatedly at minimal computational cost. The ability to simulate a large number of samples from the posterior distribution results in very accurate estimates of the inverse solution and its associated uncertainty. Combined with high accuracy of the gPC-based forward solver, the new algorithm can provide great efficiency in practical applications. A rigorous error analysis of the algorithm is conducted, where we establish convergence of the approximate posterior to the true posterior and obtain an estimate of the convergence rate. It is proved that fast (exponential) convergence of the gPC forward solution yields similarly fast (exponential) convergence of the posterior. The numerical strategy and the predicted convergence rates are then demonstrated on nonlinear inverse problems ofvarying smoothness and dimension.

BibTeX Citation

@article{ art_9,
  title = { A stochastic collocation approach to Bayesian inference in inverse problems },
  editor = {  },
  author = { Y. M. Marzouk and D. Xiu },
  journal = { Communications in Computational Physics },
  doi = { DOI:prism/16 },
  year = { 2009 },
  volume = { 6 },
  number = { 1 },
  month = { 0 },
  pages = { 826-847 },
}

Available Links

• http://dx.doi.org/10.1166/jctn.2009.1285

Abstract

Stochastic reaction networks modeled as jump Markov processes serve as the main mathematical representation of biochemical phenomena in cells, particularly when the relevant molecule count is low, causing deterministic macroscale chemical reaction models to fail. Further, as there is mainly empirical knowledge about the rate parameters, parametric uncertainty analysis becomes very important. The conventional predictability tools for deterministic systems do not readily generalize to the stochastic setting. We use spectral polynomial chaos expansions to represent stochastic processes. Bayesian inference techniques with Markov chain Monte Carlo are used to find the best spectral representation of the system state, taking into account not only intrinsic stochastic noise but also parametric uncertainties. A likelihood-based adaptive domain decomposition is introduced and applied, in particular, for the cases when the parameter range includes deterministic bifurcations. We show that the adaptive multidomain polynomial chaos representation captures the correct system behavior for a benchmark bistable Schlögl model for a wide range of parameter variations.

BibTeX Citation

@article{ art_7,
  title = { Bayesian inference of spectral expansions for predictability assessment in stochastic reaction networks },
  keywords = { uncertainty quantification; bayesian inference; polynomial chaos; stochastic reaction networks; domain decomposition; predictability },
  editor = {  },
  author = { K. Sargsyan and B. J. Debusschere and H. N. Najm and Y. M. Marzouk },
  journal = { Journal of Computational and Theoretical Nanoscience },
  doi = { 10.1166/jctn.2009.1285 },
  year = { 2009 },
  volume = { 6 },
  number = { 10 },
  month = { 0 },
  pages = { 2283-2297 },
}

Available Links

• http://dx.doi.org/10.1017/S0022112006004411

Abstract

Transverse jets arise in many applications, including propulsion, effluent dispersion, oil field flows, and V/STOL aerodynamics. This study seeks a fundamental, mechanistic understanding of the structure and evolution of vorticity in the transverse jet. We develop a high-resolution three-dimensional vortex simulation of the transverse jet at large Reynolds number and consider jet-to-crossflow velocity ratios r ranging from 5 to 10. A new formulation of vorticity-flux boundary conditions accounts for the interaction of channel wall vorticity with the jet flow immediately around the orifice. We demonstrate that the nascent jet shear layer contains not only azimuthal vorticity generated in the jet pipe, but wall-normal and azimuthal perturbations resulting from the jet--crossflow interaction. This formulation also yields analytical expressions for vortex lines in the near field as a function of $r$. Transformation of the cylindrical shear layer emanating from the orifice begins with axial elongation of its lee side to form sections of counter-rotating vorticity aligned with the jet trajectory. Periodic roll-up of the shear layer accompanies this deformation, creating complementary vortex arcs on the lee and windward sides of the jet. Counter-rotating vorticity then drives lee-side roll-ups in the windward direction, along the normal to the jet trajectory. Azimuthal vortex arcs of alternating sign thus approach each other on the windward boundary of the jet. Accordingly, initially planar material rings on the shear layer fold completely and assume an interlocking structure that persists for several diameters above the jet exit. Though the near field of the jet is dominated by deformation and periodic roll-up of the shear layer, the resulting counter-rotating vorticity is a pronounced feature of the mean field; in turn, the mean counter-rotation exerts a substantial influence on the deformation of the shear layer. Following the pronounced bending of the trajectory into the crossflow, we observe a sudden breakdown of near-field vortical structures into a dense distribution of smaller scales. Spatial filtering of this region reveals the persistence of counter-rotating streamwise vorticity initiated in the near field.

BibTeX Citation

@article{ art_5,
  title = { Vorticity structure and evolution in a transverse jet },
  editor = {  },
  author = { Y. M. Marzouk and A. F. Ghoniem },
  journal = { Journal of Fluid Mechanics },
  doi = { 10.1017/S0022112006004411 },
  year = { 2007 },
  volume = { 575 },
  number = { 1 },
  month = { 0 },
  pages = { 267-305 },
}

Available Links

• http://dx.doi.org/10.2514/2.2018

Abstract

Combustion in turbulent flows may take the form of a thin flame wrapped around vortical structures. For this regime, the flame embedding approach seeks to decouple computations of the outer nonreacting flow and the combustion zone by discretizing the flame surface into a number of elemental flames, each incorporating the local impact of unsteady flow-flame interaction. An unsteady strained laminar flame solver, based on a boundary-layer approximation of combustion in a time-dependent stagnation-point potential flow, is proposed as an elemental flame model. To validate the concept, two-dimensional simulations of premixed flame-vortex interactions are performed for a matrix of vortex strengths and length scales, and a section of the flame is selected for comparison with the flame embedding model results. Results show that using the flame leading-edge strain rate gives reasonable agreement in the cases of low strain rate and weak strain rate gradient within the flame structure. This agreement deteriorates substantially when both are high. We propose two different schemes, both based on averaging the strain rate across the flame structure, and demonstrate that agreement between the one-dimensional model and the two-dimensional simulation greatly improves when the actual strain rate at the reaction zone of the one-dimensional flame is made to match that of the two-dimensional flame.

BibTeX Citation

@article{ art_3,
  keywords = { Two dimensional model; One dimensional model; Premixed flame; Finite difference method; Numerical simulation; Flame structure; Turbulent flame; Combustion },
  title = { Toward a flame embedding model for turbulent combustion simulation },
  author = { Y. M. Marzouk and A. F. Ghoniem and H. N. Najm },
  journal = { AIAA Journal },
  volume = { 41 },
  number = { 4 },
  pages = { 641-652 },
  month = { 0 },
  year = { 2003 },
  editor = {  },
}

Available Links

• http://dx.doi.org/10.1007/s003480050247

Abstract

Autocorrelation of a double-exposed image, unlike cross-correlation between two images, produces a correlation function that is symmetric about the origin. Thus, while it is possible to calculate the speed and direction of tracer particles in a particle image velocimetry (PIV) image using autocorrelation, it is impossible to tell whether the velocity is in the positive or negative direction. This ambiguity can be resolved by spatially shifting one exposure relative to the next or labeling exposures with color or polarization, but the complexity and limitations of these methods can be prohibitive. It is, however, possible to resolve the sign of the velocity from a triple-exposed image using unequal time intervals between exposures. Triple-exposed images, like double-exposed images, correlate symmetrically about zero. The directional ambiguity, however, can be resolved by calculating the probability that the three exposures occur in a specific temporal order; that is, by assuming that the correlation has a specific sign and testing to see if the assumption is correct. Traditional spectral and statistical correlation techniques are unable to accomplish this. Presented herein is a computationally efficient asymmetric correlation function that is able to differentiate the temporal order of triple exposed images. Included is a discussion of the limitations of this function and of difficulties in experimental implementation.

BibTeX Citation

@article{ art_1,
  title = { Asymmetric autocorrelation function to resolve directional ambiguity in PIV images },
  editor = {  },
  author = { Y. M. Marzouk and D. P. Hart },
  journal = { Experiments in Fluids },
  doi = { 10.1007/s003480050247 },
  year = { 1998 },
  volume = { 25 },
  number = { 5-6 },
  month = { 0 },
  pages = { 401-408 },
}