MIT Uncertainty Quantification Group

@article{li-sharpdetectionlowdimensional-2025,
  author     = "Li, M. T. C. and Cui, T. and Li, F. and Marzouk, Y. M. and Zahm, O.",
  journal    = "Information and Inference: A Journal of the IMA",
  number     = "3",
  publisher  = "Oxford University Press",
  title      = "Sharp detection of low-dimensional structure in probability measures via dimensional logarithmic Sobolev inequalities",
  volume     = "14",
  year       = "2025"
}

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

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.

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

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.

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

Data generation remains a bottleneck in training surrogate models to predict molecular properties. We demonstrate that multitask Gaussian process regression overcomes this limitation by leveraging both expensive and cheap data sources. In particular, we consider training sets constructed from coupled-cluster (CC) and density functional theory (DFT) data. We report that multitask surrogates can predict at CC-level accuracy with a reduction in data generation cost by over an order of magnitude. Of note, our approach allows the training set to include DFT data generated by a heterogeneous mix of exchange–correlation functionals without imposing any artificial hierarchy on functional accuracy. More generally, the multitask framework can accommodate a wider range of training set structures—including the full disparity between the different levels of fidelity—than existing kernel approaches based on Δ-learning although we show that the accuracy of the two approaches can be similar. Consequently, multitask regression can be a tool for reducing data generation costs even further by opportunistically exploiting existing data sources.

@article{fisher-multitask-2024,
  author     = "Fisher, K. E. and Herbst, M. F. and Marzouk, Y. M.",
  journal    = "The Journal of Chemical Physics",
  month      = "7",
  number     = "014114",
  title      = "Multitask methods for predicting molecular properties from heterogeneous data",
  volume     = "161",
  year       = "2024",
  doi        = "10.1063/5.0201681",
  abstract   = "Data generation remains a bottleneck in training surrogate models to predict molecular properties. We demonstrate that multitask Gaussian process regression overcomes this limitation by leveraging both expensive and cheap data sources. In particular, we consider training sets constructed from coupled-cluster (CC) and density functional theory (DFT) data. We report that multitask surrogates can predict at CC-level accuracy with a reduction in data generation cost by over an order of magnitude. Of note, our approach allows the training set to include DFT data generated by a heterogeneous mix of exchange–correlation functionals without imposing any artificial hierarchy on functional accuracy. More generally, the multitask framework can accommodate a wider range of training set structures—including the full disparity between the different levels of fidelity—than existing kernel approaches based on Δ-learning although we show that the accuracy of the two approaches can be similar. Consequently, multitask regression can be a tool for reducing data generation costs even further by opportunistically exploiting existing data sources.",
  keywords   = "Quantum chemistry, Coupled-cluster methods, Density functional theory, Exchange correlation functionals, Machine learning, Molecular properties, Covariance and correlation, Gaussian processes, Nonparametric statistics"
}

We propose a novel diffusion map particle system (DMPS) for generative modeling, based on diffusion maps and Laplacian-adjusted Wasserstein gradient descent (LAWGD). Diffusion maps are used to approximate the generator of the Langevin diffusion process from samples, and hence to learn the underlying data-generating manifold. On the other hand, LAWGD enables efficient sampling from the target distribution given a suitable choice of kernel, which we construct here via a spectral approximation of the generator, computed with diffusion maps. Numerical experiments show that our method outperforms others on synthetic datasets, including examples with manifold structure.

@article{li-diffusionparticles-2023,
  author     = "Li, F. Y. and Marzouk, Y. M.",
  booktitle  = "ICML 2023 Workshop on Structured Probabilistic Inference & Generative Modeling",
  journal    = "Foundations of Data Science",
  month      = "9",
  number     = "3",
  pages      = "814--837",
  title      = "Diffusion map particle systems for generative modeling",
  volume     = "7",
  year       = "2025",
  doi        = "10.3934/fods.2024054",
  abstract   = "We propose a novel diffusion map particle system (DMPS) for generative modeling, based on diffusion maps and Laplacian-adjusted Wasserstein gradient descent (LAWGD). Diffusion maps are used to approximate the generator of the Langevin diffusion process from samples, and hence to learn the underlying data-generating manifold. On the other hand, LAWGD enables efficient sampling from the target distribution given a suitable choice of kernel, which we construct here via a spectral approximation of the generator, computed with diffusion maps. Numerical experiments show that our method outperforms others on synthetic datasets, including examples with manifold structure.",
  keywords   = "diffusion maps, kernel methods, gradient flows, generative modeling, sampling"
}

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.

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

The joint retrieval of surface reflectances and atmospheric parameters in visible/short-wave infrared (VSWIR) imaging spectroscopy is a computationally challenging high-dimensional problem. Using NASA’s Surface Biology and Geology mission (SBG) as the motivational context, the uncertainty associated with the retrievals is crucial for further application of the retrieved results for environmental applications. Although Markov chain Monte Carlo (MCMC) is a Bayesian method ideal for uncertainty quantification (UQ), the full-dimensional implementation of MCMC for the retrieval is computationally intractable. In this work, we developed a block Metropolis MCMC algorithm for the high-dimensional VSWIR surface reflectance retrieval that leverages conditional linearity structure in the forward radiative transfer model to enable tractable fully Bayesian computation. We use the posterior distribution from this MCMC algorithm to assess the limitations of optimal estimation (OE), the state-of-the-art Bayesian algorithm in operational retrievals which is more computationally efficient but uses a Gaussian approximation to characterize the posterior. Analyzing the differences in the posterior computed by each method, the MCMC algorithm was shown to give more physically sensible results and reveals the non-Gaussian structure of the posterior, specifically in the atmospheric aerosol optical depth (AOD) parameter and the low-wavelength surface reflectances.

@article{leung-spectroscopy-2023,
  author     = "Leung, K. and Thompson, D. and Susiluoto, J. and Jagalur, J. and Braverman, A. and Marzouk, Y. M.",
  journal    = "IEEE Transactions on Geoscience and Remote Sensing",
  month      = "6",
  pages      = "1-12",
  title      = "Evaluating the accuracy of Gaussian approximations in VSWIR imaging spectroscopy retrievals",
  volume     = "62",
  year       = "2024",
  doi        = "10.1109/TGRS.2024.3411916",
  abstract   = "The joint retrieval of surface reflectances and atmospheric parameters in visible/short-wave infrared (VSWIR) imaging spectroscopy is a computationally challenging high-dimensional problem. Using NASA’s Surface Biology and Geology mission (SBG) as the motivational context, the uncertainty associated with the retrievals is crucial for further application of the retrieved results for environmental applications. Although Markov chain Monte Carlo (MCMC) is a Bayesian method ideal for uncertainty quantification (UQ), the full-dimensional implementation of MCMC for the retrieval is computationally intractable. In this work, we developed a block Metropolis MCMC algorithm for the high-dimensional VSWIR surface reflectance retrieval that leverages conditional linearity structure in the forward radiative transfer model to enable tractable fully Bayesian computation. We use the posterior distribution from this MCMC algorithm to assess the limitations of optimal estimation (OE), the state-of-the-art Bayesian algorithm in operational retrievals which is more computationally efficient but uses a Gaussian approximation to characterize the posterior. Analyzing the differences in the posterior computed by each method, the MCMC algorithm was shown to give more physically sensible results and reveals the non-Gaussian structure of the posterior, specifically in the atmospheric aerosol optical depth (AOD) parameter and the low-wavelength surface reflectances.",
  keywords   = "mcmc"
}

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.

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

We consider the problem of sampling from an unknown distribution for which only a sufficiently large number of training samples are available. Such settings have recently drawn considerable interest in the context of generative modelling and Bayesian inference. In this paper, we propose a generative model combining Schrödinger bridges and Langevin dynamics. Schrödinger bridges over an appropriate reversible reference process are used to approximate the conditional transition probability from the available training samples, which is then implemented in a discrete-time reversible Langevin sampler to generate new samples. By setting the kernel bandwidth in the reference process to match the time step size used in the unadjusted Langevin algorithm, our method effectively circumvents any stability issues typically associated with the time-stepping of stiff stochastic differential equations. Moreover, we introduce a novel split-step scheme, ensuring that the generated samples remain within the convex hull of the training samples. Our framework can be naturally extended to generate conditional samples and to Bayesian inference problems. We demonstrate the performance of our proposed scheme through experiments on synthetic datasets with increasing dimensions and on a stochastic subgrid-scale parametrization conditional sampling problem.

@article{li-schroedinger-2024,
  author     = "Gottwald, G. and Li, F. Y. and Marzouk, Y. M. and Reich, S.",
  journal    = "Philosophical Transactions of the Royal Society A",
  month      = "6",
  number     = "2299",
  pages      = "20240332",
  title      = "Stable generative modeling using Schrödinger bridges",
  volume     = "383",
  year       = "2025",
  doi        = "https://doi.org/10.48550/arXiv.2401.04372",
  abstract   = "We consider the problem of sampling from an unknown distribution for which only a sufficiently large number of training samples are available. Such settings have recently drawn considerable interest in the context of generative modelling and Bayesian inference. In this paper, we propose a generative model combining Schrödinger bridges and Langevin dynamics. Schrödinger bridges over an appropriate reversible reference process are used to approximate the conditional transition probability from the available training samples, which is then implemented in a discrete-time reversible Langevin sampler to generate new samples. By setting the kernel bandwidth in the reference process to match the time step size used in the unadjusted Langevin algorithm, our method effectively circumvents any stability issues typically associated with the time-stepping of stiff stochastic differential equations. Moreover, we introduce a novel split-step scheme, ensuring that the generated samples remain within the convex hull of the training samples. Our framework can be naturally extended to generate conditional samples and to Bayesian inference problems. We demonstrate the performance of our proposed scheme through experiments on synthetic datasets with increasing dimensions and on a stochastic subgrid-scale parametrization conditional sampling problem.",
  keywords   = "Langevin dynamics, Schrödinger bridges, generative modeling, Bayesian inference"
}

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.

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

This article presents a general approximation-theoretic framework to analyze measure transport algorithms for probabilistic modeling. A primary motivating application for such algorithms is sampling—a central task in statistical inference and generative modeling. We provide a priori error estimates in the continuum limit, i.e., when the measures (or their densities) are given, but when the transport map is discretized or approximated using a finite-dimensional function space. Our analysis relies on the regularity theory of transport maps and on classical approximation theory for high-dimensional functions. A third element of our analysis, which is of independent interest, is the development of new stability estimates that relate the distance between two maps to the distance (or divergence) between the pushforward measures they define. We present a series of applications of our framework, where quantitative convergence rates are obtained for practical problems using Wasserstein metrics, maximum mean discrepancy, and Kullback–Leibler divergence. Specialized rates for approximations of the popular triangular Knöthe–Rosenblatt maps are obtained, followed by numerical experiments that demonstrate and extend our theory.

@article{baptista-transportapproxtheory-2024,
  author     = "Baptista, R. and Hosseini, B. and Kovachki, N. B. and Marzouk, Y. M. and Sagiv, A.",
  journal    = "Mathematics of Computation",
  title      = "An approximation theory framework for measure-transport sampling algorithms",
  year       = "2025",
  doi        = "10.1090/mcom/4013",
  abstract   = "This article presents a general approximation-theoretic framework to analyze measure transport algorithms for probabilistic modeling. A primary motivating application for such algorithms is sampling—a central task in statistical inference and generative modeling. We provide a priori error estimates in the continuum limit, i.e., when the measures (or their densities) are given, but when the transport map is discretized or approximated using a finite-dimensional function space. Our analysis relies on the regularity theory of transport maps and on classical approximation theory for high-dimensional functions. A third element of our analysis, which is of independent interest, is the development of new stability estimates that relate the distance between two maps to the distance (or divergence) between the pushforward measures they define. We present a series of applications of our framework, where quantitative convergence rates are obtained for practical problems using Wasserstein metrics, maximum mean discrepancy, and Kullback–Leibler divergence. Specialized rates for approximations of the popular triangular Knöthe–Rosenblatt maps are obtained, followed by numerical experiments that demonstrate and extend our theory.",
  keywords   = "Transport map, generative models, stability analysis, approximation theory"
}

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.

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

We present a novel framework for conditional sampling of probability measures, using block triangular transport maps. We develop the theoretical foundations of block triangular transport in a Banach space setting, establishing general conditions under which conditional sampling can be achieved and drawing connections between monotone block triangular maps and optimal transport. Based on this theory, we then introduce a computational approach, called monotone generative adversarial networks (M-GANs), to learn suitable block triangular maps. Our algorithm uses only samples from the underlying joint probability measure and is hence likelihood-free. Numerical experiments with M-GAN demonstrate accurate sampling of conditional measures in synthetic examples, Bayesian inverse problems involving ordinary and partial differential equations, and probabilistic image inpainting.

@article{baptista-monotonegans-2024,
  author     = "Baptista, R. and Hosseini, B. and Kovachki, N. B. and Marzouk, Y. M.",
  journal    = "SIAM/ASA Journal on Uncertainty Quantification",
  month      = "8",
  number     = "3",
  pages      = "868-900",
  title      = "Conditional Sampling with Monotone GANs: From Generative Models to Likelihood-Free Inference",
  volume     = "12",
  year       = "2024",
  doi        = "10.1137/23M1581546",
  abstract   = "We present a novel framework for conditional sampling of probability measures, using block triangular transport maps. We develop the theoretical foundations of block triangular transport in a Banach space setting, establishing general conditions under which conditional sampling can be achieved and drawing connections between monotone block triangular maps and optimal transport. Based on this theory, we then introduce a computational approach, called monotone generative adversarial networks (M-GANs), to learn suitable block triangular maps. Our algorithm uses only samples from the underlying joint probability measure and is hence likelihood-free. Numerical experiments with M-GAN demonstrate accurate sampling of conditional measures in synthetic examples, Bayesian inverse problems involving ordinary and partial differential equations, and probabilistic image inpainting.",
  keywords   = "measure transport;conditional simulation;likelihood-free inference;optimal transport;GANs;normalizing flows"
}

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.

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

Transportation of measure provides a versatile approach for modeling complex probability distributions, with applications in density estimation, Bayesian inference, generative modeling, and beyond. Monotone triangular transport maps—approximations of the Knothe–Rosenblatt (KR) rearrangement—are a canonical choice for these tasks. Yet the representation and parameterization of such maps have a significant impact on their generality and expressiveness, and on properties of the optimization problem that arises in learning a map from data (e.g., via maximum likelihood estimation). We present a general framework for representing monotone triangular maps via invertible transformations of smooth functions. We establish conditions on the transformation such that the associated infinite-dimensional minimization problem has no spurious local minima, i.e., all local minima are global minima; and we show for target distributions satisfying certain tail conditions that the unique global minimizer corresponds to the KR map. Given a sample from the target, we then propose an adaptive algorithm that estimates a sparse semi-parametric approximation of the underlying KR map. We demonstrate how this framework can be applied to joint and conditional density estimation, likelihood-free inference, and structure learning of directed graphical models, with stable generalization performance across a range of sample sizes.

@article{baptista-atm-2023,
  author     = "Baptista, R. and Marzouk, Y. M. and Zahm, O.",
  journal    = "Foundations of Computational Mathematics",
  month      = "11",
  title      = "On the representation and learning of monotone triangular transport maps",
  year       = "2023",
  doi        = "10.1007/s10208-023-09630-x",
  abstract   = "Transportation of measure provides a versatile approach for modeling complex probability distributions, with applications in density estimation, Bayesian inference, generative modeling, and beyond. Monotone triangular transport maps—approximations of the Knothe–Rosenblatt (KR) rearrangement—are a canonical choice for these tasks. Yet the representation and parameterization of such maps have a significant impact on their generality and expressiveness, and on properties of the optimization problem that arises in learning a map from data (e.g., via maximum likelihood estimation). We present a general framework for representing monotone triangular maps via invertible transformations of smooth functions. We establish conditions on the transformation such that the associated infinite-dimensional minimization problem has no spurious local minima, i.e., all local minima are global minima; and we show for target distributions satisfying certain tail conditions that the unique global minimizer corresponds to the KR map. Given a sample from the target, we then propose an adaptive algorithm that estimates a sparse semi-parametric approximation of the underlying KR map. We demonstrate how this framework can be applied to joint and conditional density estimation, likelihood-free inference, and structure learning of directed graphical models, with stable generalization performance across a range of sample sizes.",
  keywords   = "Transportation of measure, Knothe–Rosenblatt rearrangement, normalizing flows, monotone functions, infinite-dimensional optimization, adaptive approximation, multivariate polynomials, wavelets, density estimation."
}

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.

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

Smoothing is a specialized form of Bayesian inference for state-space models that characterizes the posterior distribution of a collection of states given an associated sequence of observations. Our companion manuscript proposes a general framework for transport-based ensemble smoothing, which includes linear Kalman-type smoothers as special cases. Here, we build on this foundation to realize and demonstrate nonlinear backward ensemble transport smoothers. We discuss parameterization and regularization of the associated transport maps, and then examine the performance of these smoothers for nonlinear and chaotic dynamical systems that exhibit non-Gaussian behavior. In these settings, our nonlinear transport smoothers yield lower estimation error than conventional linear smoothers and state-of-the-art iterative ensemble Kalman smoothers, for comparable numbers of model evaluations.

@article{ramgraber-smoothingparttwo-2022,
  author     = "Ramgraber, M. and Baptista, R. and McLaughlin, D. and Marzouk, Y. M.",
  journal    = "Journal of Computational Physics: X",
  month      = "11",
  number     = "100133",
  title      = "Ensemble transport smoothing. Part 2: nonlinear updates",
  volume     = "17",
  year       = "2023",
  doi        = "10.1016/j.jcpx.2023.100133",
  abstract   = "Smoothing is a specialized form of Bayesian inference for state-space models that characterizes the posterior distribution of a collection of states given an associated sequence of observations. Our companion manuscript proposes a general framework for transport-based ensemble smoothing, which includes linear Kalman-type smoothers as special cases. Here, we build on this foundation to realize and demonstrate nonlinear backward ensemble transport smoothers. We discuss parameterization and regularization of the associated transport maps, and then examine the performance of these smoothers for nonlinear and chaotic dynamical systems that exhibit non-Gaussian behavior. In these settings, our nonlinear transport smoothers yield lower estimation error than conventional linear smoothers and state-of-the-art iterative ensemble Kalman smoothers, for comparable numbers of model evaluations.",
  keywords   = "Data assimilation, smoothing, ensemble methods, triangular transport."
}

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.

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

@article{parno-mpart-2022,
  author     = "Parno, M. and Rubio, P.-B. and Sharp, D. and Brennan, M. and Baptista, R. and Bonart, H. and Marzouk, Y. M.",
  journal    = "Journal of Open Source Software",
  number     = "80",
  pages      = "4843",
  title      = "MParT: Monotone Parameterization Toolkit",
  volume     = "7",
  year       = "2022"
}

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.

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

We introduce a method for the nonlinear dimension reduction of a high-dimensional function $u: \mathbb{R}^d \to \mathbb{R}$, $d \gg 1$. Our objective is to identify a nonlinear feature map $g: \mathbb{R}^d \to \mathbb{R}^m$, with a prescribed intermediate dimension $m \ll d$, so that $u$ can be well approximated by $f \circ g$ for some profile function $f: \mathbb{R}^m \to \mathbb{R}$. We propose to build the feature map by aligning the Jacobian $\nabla g$ with the gradient $\nabla u$, and we theoretically analyze the properties of the resulting $g$. Once $g$ is built, we construct $f$ by solving a gradient-enhanced least squares problem. Our practical algorithm makes use of a sample $\{ x(i), u(x(i)), \nabla u(x(i)) \}_{i=1}^N$ and builds both $g$ and $f$ on adaptive downward-closed polynomial spaces, using cross validation to avoid overfitting. We numerically evaluate the performance of our algorithm across different benchmarks, and explore the impact of the intermediate dimension $m$. We show that building a nonlinear feature map $g$ can permit more accurate approximation of $u$ than a linear $g$, for the same input data set.

@article{bigoni-nonlineardimred-2021,
  author     = "Bigoni, D. and Marzouk, Y. M. and Prieur, C. and Zahm, O.",
  journal    = "Information and Inference: A Journal of the IMA",
  title      = "Nonlinear dimension reduction for surrogate modeling using gradient information",
  volume     = "in press",
  year       = "2022",
  doi        = "10.1093/imaiai/iaac006",
  abstract   = "We introduce a method for the nonlinear dimension reduction of a high-dimensional function $u: \mathbb{R}^d \to \mathbb{R}$, $d \gg 1$. Our objective is to identify a nonlinear feature map $g: \mathbb{R}^d  \to \mathbb{R}^m$, with a prescribed intermediate dimension $m \ll d$, so that $u$ can be well approximated by $f \circ g$ for some profile function $f: \mathbb{R}^m \to \mathbb{R}$. We propose to build the feature map by aligning the Jacobian $\nabla g$ with the gradient $\nabla u$, and we theoretically analyze the properties of the resulting $g$. Once $g$ is built, we construct $f$ by solving a gradient-enhanced least squares problem. Our practical algorithm makes use of a sample $\{ x(i), u(x(i)), \nabla u(x(i)) \}_{i=1}^N$ and builds both $g$ and $f$ on adaptive downward-closed polynomial spaces, using cross validation to avoid overfitting. We numerically evaluate the performance of our algorithm across different benchmarks, and explore the impact of the intermediate dimension $m$. We show that building a nonlinear feature map $g$ can permit more accurate approximation of $u$ than a linear $g$, for the same input data set."
}

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.

@article{art_79,
  author     = "Spantini, A. and Baptista, R. and Marzouk, Y. M.",
  journal    = "SIAM Review",
  number     = "4",
  pages      = "921–953",
  title      = "Coupling techniques for nonlinear ensemble filtering",
  volume     = "64",
  year       = "2022",
  doi        = "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.",
  keywords   = "nonlinear filtering, state-space models, couplings, transport maps, ensemble Kalman filter, graphical models, localization, approximate Bayesian computation"
}

We propose a regularization method for ensemble Kalman filtering (EnKF) with elliptic observation operators. Commonly used EnKF regularization methods suppress state correlations at long distances. For observations described by elliptic partial differential equations, such as the pressure Poisson equation (PPE) in incompressible fluid flows, distance localization should be used cautiously, as we cannot disentangle slowly decaying physical interactions from spurious long-range correlations. This is particularly true for the PPE, in which distant vortex elements couple nonlinearly to induce pressure. Instead, these inverse problems have a low effective dimension: low-dimensional projections of the observations strongly inform a low-dimensional subspace of the state space. We derive a low-rank factorization of the Kalman gain based on the spectrum of the Jacobian of the observation operator. The identified eigenvectors generalize the source and target modes of the multipole expansion, independently of the underlying spatial distribution of the problem. Given rapid spectral decay, inference can be performed in the low-dimensional subspace spanned by the dominant eigenvectors. This low-rank EnKF is assessed on dynamical systems with Poisson observation operators, where we seek to estimate the positions and strengths of point singularities over time from potential or pressure observations. We also comment on the broader applicability of this approach to elliptic inverse problems outside the context of filtering.

@article{leprovost-lrenkf-2022,
  author     = "Le Provost, M. and Baptista, R. and Marzouk, Y. M. and Eldredge, J.",
  journal    = "Proceedings of the Royal Society A",
  number     = "2266",
  title      = "A low-rank ensemble Kalman filter for elliptic observations",
  volume     = "478",
  year       = "2022",
  doi        = "10.1098/rspa.2022.0182",
  abstract   = "We propose a regularization method for ensemble Kalman filtering (EnKF) with elliptic observation operators. Commonly used EnKF regularization methods suppress state correlations at long distances. For observations described by elliptic partial differential equations, such as the pressure Poisson equation (PPE) in incompressible fluid flows, distance localization should be used cautiously, as we cannot disentangle slowly decaying physical interactions from spurious long-range correlations. This is particularly true for the PPE, in which distant vortex elements couple nonlinearly to induce pressure. Instead, these inverse problems have a low effective dimension: low-dimensional projections of the observations strongly inform a low-dimensional subspace of the state space. We derive a low-rank factorization of the Kalman gain based on the spectrum of the Jacobian of the observation operator. The identified eigenvectors generalize the source and target modes of the multipole expansion, independently of the underlying spatial distribution of the problem. Given rapid spectral decay, inference can be performed in the low-dimensional subspace spanned by the dominant eigenvectors. This low-rank EnKF is assessed on dynamical systems with Poisson observation operators, where we seek to estimate the positions and strengths of point singularities over time from potential or pressure observations. We also comment on the broader applicability of this approach to elliptic inverse problems outside the context of filtering.",
  keywords   = "data assimilation, ensemble Kalman filter, elliptic inverse problems, incompressible fluid mechanics, observation-informed dimension reduction"
}

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.

@article{litvinenko-tensordiv-2022,
  author     = "Litvinenko, A. and Marzouk, Y. M. and Matthies, H. G. and Scavino, M. and Spantini, A.",
  journal    = "Numerical Linear Algebra with Applications",
  title      = "Computing f-divergences and distances of high-dimensional probability density functions",
  volume     = "in press",
  year       = "2022",
  doi        = "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."
}

We develop several inference methods to estimate the position of dislocations from images generated using dark-field X-ray microscopy (DFXM)—achieving superresolution accuracy and principled uncertainty quantification. Using the framework of Bayesian inference, we incorporate models of the DFXM contrast mechanism and detector measurement noise, along with initial position estimates, into a statistical model coupling DFXM images with the dislocation position of interest. We motivate several position estimation and uncertainty quantification algorithms based on this model. We then demonstrate the accuracy of our primary estimation algorithm on synthetic realistic DFXM images of edge dislocations in single-crystal aluminum. We conclude with a discussion of our methods’ impact on future dislocation studies and possible future research avenues.

@article{brennan-dfxm-2022,
  author     = "Brennan, M. and Howard, M. and Marzouk, Y. M. and Dresselhaus-Marais, L.",
  journal    = "Journal of Materials Science",
  pages      = "14890–14904",
  title      = "Analytical methods for superresolution dislocation identification in dark-field X-ray microscopy",
  volume     = "57",
  year       = "2022",
  doi        = "10.1007/s10853-022-07465-5",
  abstract   = "We develop several inference methods to estimate the position of dislocations from images generated using dark-field X-ray microscopy (DFXM)—achieving superresolution accuracy and principled uncertainty quantification. Using the framework of Bayesian inference, we incorporate models of the DFXM contrast mechanism and detector measurement noise, along with initial position estimates, into a statistical model coupling DFXM images with the dislocation position of interest. We motivate several position estimation and uncertainty quantification algorithms based on this model. We then demonstrate the accuracy of our primary estimation algorithm on synthetic realistic DFXM images of edge dislocations in single-crystal aluminum. We conclude with a discussion of our methods’ impact on future dislocation studies and possible future research avenues.",
  keywords   = "dislocation, Bayesian inference, dark-field X-ray microscopy, metals"
}

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.

@article{zhang-girr-2021,
  author     = "Zhang, B. J. and Marzouk, Y. M. and Spiliopoulos, K.",
  journal    = "Statistics and Computing",
  number     = "78",
  title      = "Geometry-informed irreversible perturbations for accelerated convergence of Langevin dynamics",
  volume     = "32",
  year       = "2022",
  doi        = "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."
}

We propose a dimension reduction technique for Bayesian inverse problems with nonlinear forward operators, non-Gaussian priors, and non-Gaussian observation noise. The likelihood function is approximated by a ridge function, i.e., a map which depends non-trivially only on a few linear combinations of the parameters. We build this ridge approximation by minimizing an upper bound on the Kullback–Leibler divergence between the posterior distribution and its approximation. This bound, obtained via logarithmic Sobolev inequalities, allows one to certify the error of the posterior approximation. Computing the bound requires computing the second moment matrix of the gradient of the log likelihood function. In practice, a sample-based approximation of the upper bound is then required. We provide an analysis that enables control of the posterior approximation error due to this sampling. Numerical and theoretical comparisons with existing methods illustrate the benefits of the proposed methodology.

@article{art_73,
  author     = "Zahm, O. and Cui, T. and Law, K. J. H. and Spantini, A. and Marzouk, Y. M.",
  journal    = "Mathematics of Computation",
  pages      = "1789–1835",
  title      = "Certified dimension reduction in nonlinear Bayesian inverse problems",
  volume     = "91",
  year       = "2022",
  doi        = "10.1090/mcom/3737",
  abstract   = "We propose a dimension reduction technique for Bayesian inverse problems with nonlinear forward operators, non-Gaussian priors, and non-Gaussian observation noise. The likelihood function is approximated by a ridge function, i.e., a map which depends non-trivially only on a few linear combinations of the parameters. We build this ridge approximation by minimizing an upper bound on the Kullback--Leibler divergence between the posterior distribution and its approximation. This bound, obtained via logarithmic Sobolev inequalities, allows one to certify the error of the posterior approximation. Computing the bound requires computing the second moment matrix of the gradient of the log likelihood function. In practice, a sample-based approximation of the upper bound is then required. We provide an analysis that enables control of the posterior approximation error due to this sampling. Numerical and theoretical comparisons with existing methods illustrate the benefits of the proposed methodology.",
  keywords   = "dimension reduction, nonlinear Bayesian inverse problem, logarithmic Sobolev inequality, certified error bound, non-asymptotic analysis"
}

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.

@article{zech-infdim-2021,
  author     = "Zech, J. and Marzouk, Y. M.",
  journal    = "Constructive Approximation",
  pages      = "987–1036",
  title      = "Sparse approximation of triangular transports. Part II: the infinite dimensional case",
  volume     = "55",
  year       = "2022",
  doi        = "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."
}

For two probability measures $\rho$ and $\pi$ with analytic densities on the $d$-dimensional cube $[−1,1]^d$, we investigate the approximation of the unique triangular monotone Knothe–Rosenblatt transport $T: [−1,1]^d \to [−1,1]^d$, such that the pushforward $T_\sharp \rho$ equals $\pi$. It is shown that for $d \in \mathbb{N}$ there exist approximations $\tilde{T}$ of $T$, based on either sparse polynomial expansions or deep ReLU neural networks, such that the distance between $\tilde{T}_\sharp \rho$ and $\pi$ decreases exponentially. More precisely, we prove error bounds of the type $\exp(−\beta N^{1/d})$ (or $\exp(−\beta N^{1/d})$) for neural networks), where $N$ refers to the dimension of the ansatz space (or the size of the network) containing $\tilde{T}$; the notion of distance comprises the Hellinger distance, the total variation distance, the Wasserstein distance and the Kullback–Leibler divergence. Our construction guarantees $\tilde{T}$ to be a monotone triangular bijective transport on the hypercube $[−1,1]^d$. Analogous results hold for the inverse transport $S=T^{−1}$. The proofs are constructive, and we give an explicit a priori description of the ansatz space, which can be used for numerical implementations.

@article{zech-triangular-2020,
  author     = "Zech, J. and Marzouk, Y. M.",
  journal    = "Constructive Approximation",
  pages      = "919–986",
  title      = "Sparse approximation of triangular transports. Part I: the finite dimensional case",
  volume     = "55",
  year       = "2022",
  doi        = "10.1007/s00365-022-09569-2",
  abstract   = "For two probability measures $\rho$ and $\pi$ with analytic densities on the $d$-dimensional cube $[−1,1]^d$, we investigate the approximation of the unique triangular monotone Knothe--Rosenblatt transport $T: [−1,1]^d \to [−1,1]^d$, such that the pushforward $T_\sharp \rho$ equals $\pi$. It is shown that for $d \in \mathbb{N}$ there exist approximations $\tilde{T}$ of $T$, based on either sparse polynomial expansions or deep ReLU neural networks, such that the distance between $\tilde{T}_\sharp \rho$ and $\pi$ decreases exponentially. More precisely, we prove error bounds of the type $\exp(−\beta N^{1/d})$ (or $\exp(−\beta N^{1/d})$) for neural networks), where $N$ refers to the dimension of the ansatz space (or the size of the network) containing $\tilde{T}$; the notion of distance comprises the Hellinger distance, the total variation distance, the Wasserstein distance and the Kullback--Leibler divergence. Our construction guarantees $\tilde{T}$ to be a monotone triangular bijective transport on the hypercube $[−1,1]^d$. Analogous results hold for the inverse transport $S=T^{−1}$. The proofs are constructive, and we give an explicit a priori description of the ansatz space, which can be used for numerical implementations.",
  keywords   = "transport maps, domains of holomorphy, uncertainty quantification, sparse approximation, neural networks, sampling"
}

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.

@article{zhang-kfr-2021,
  author     = "Zhang, B. J. and Sahai, T. and Marzouk, Y. M.",
  journal    = "Journal of Computational Physics",
  pages      = "111025",
  title      = "A Koopman framework for rare event simulation in stochastic differential equations",
  volume     = "456",
  year       = "2022",
  doi        = "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."
}

Many Bayesian inference problems involve target distributions whose density functions are computationally expensive to evaluate. Replacing the target density with a local approximation based on a small number of carefully chosen density evaluations can significantly reduce the computational expense of Markov chain Monte Carlo (MCMC) sampling. Moreover, continual refinement of the local approximation can guarantee asymptotically exact sampling. We devise a new strategy for balancing the decay rate of the bias due to the approximation with that of the MCMC variance. We prove that the error of the resulting local approximation MCMC (LA-MCMC) algorithm decays at roughly the expected $1/\sqrt{T}$ rate, and we demonstrate this rate numerically. We also introduce an algorithmic parameter that guarantees convergence given very weak tail bounds, significantly strengthening previous convergence results. Finally, we apply LA-MCMC to a computationally intensive Bayesian inverse problem arising in groundwater hydrology.

@article{davis-lamcmc-2020,
  author     = "Davis, A. and Marzouk, Y. M. and Pillai, N. and Smith, A.",
  journal    = "Statistics and Computing",
  number     = "60",
  title      = "Rate-optimal refinement strategies for local approximation MCMC",
  volume     = "32",
  year       = "2022",
  doi        = "10.1007/s11222-022-10123-0",
  abstract   = "Many Bayesian inference problems involve target distributions whose density functions are computationally expensive to evaluate. Replacing the target density with a local approximation based on a small number of carefully chosen density evaluations can significantly reduce the computational expense of Markov chain Monte Carlo (MCMC) sampling. Moreover, continual refinement of the local approximation can guarantee asymptotically exact sampling. We devise a new strategy for balancing the decay rate of the bias due to the approximation with that of the MCMC variance. We prove that the error of the resulting local approximation MCMC (LA-MCMC) algorithm decays at roughly the expected $1/\sqrt{T}$ rate, and we demonstrate this rate numerically. We also introduce an algorithmic parameter that guarantees convergence given very weak tail bounds, significantly strengthening previous convergence results. Finally, we apply LA-MCMC to a computationally intensive Bayesian inverse problem arising in groundwater hydrology.",
  keywords   = "Markov chain Monte Carlo, local regression, Bayesian inference, surrogate models, sampling"
}

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.

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

We propose and analyze batch greedy heuristics for cardinality constrained maximization of non-submodular non-decreasing set functions. We consider the standard greedy paradigm, along with its distributed greedy and stochastic greedy variants. Our theoretical guarantees are characterized by the combination of submodularity and supermodularity ratios. We argue how these parameters define tight modular bounds based on incremental gains, and provide a novel reinterpretation of the classical greedy algorithm using the minorize-maximize (MM) principle. Based on that analogy, we propose a new class of methods exploiting any plausible modular bound. In the context of optimal experimental design for linear Bayesian inverse problems, we bound the submodularity and supermodularity ratios when the underlying objective is based on mutual information. We also develop novel modular bounds for the mutual information in this setting, and describe certain connections to polyhedral combinatorics. We discuss how algorithms using these modular bounds relate to established statistical notions such as leverage scores and to more recent efforts such as volume sampling. We demonstrate our theoretical findings on synthetic problems and on a real-world climate monitoring example.

@article{jagalur-batchgreedy-2020,
  author     = "Jagalur-Mohan, J. and Marzouk, Y. M.",
  journal    = "The Journal of Machine Learning Research",
  number     = "252",
  pages      = "1–62",
  title      = "Batch greedy maximization of non-submodular functions: guarantees and applications to experimental design",
  volume     = "22",
  year       = "2021",
  abstract   = "We propose and analyze batch greedy heuristics for cardinality constrained maximization of non-submodular non-decreasing set functions. We consider the standard greedy paradigm, along with its distributed greedy and stochastic greedy variants. Our theoretical guarantees are characterized by the combination of submodularity and supermodularity ratios. We argue how these parameters define tight modular bounds based on incremental gains, and provide a novel reinterpretation of the classical greedy algorithm using the minorize-maximize (MM) principle. Based on that analogy, we propose a new class of methods exploiting any plausible modular bound. In the context of optimal experimental design for linear Bayesian inverse problems, we bound the submodularity and supermodularity ratios when the underlying objective is based on mutual information. We also develop novel modular bounds for the mutual information in this setting, and describe certain connections to polyhedral combinatorics. We discuss how algorithms using these modular bounds relate to established statistical notions such as leverage scores and to more recent efforts such as volume sampling. We demonstrate our theoretical findings on synthetic problems and on a real-world climate monitoring example.",
  keywords   = "greedy methods, submodularity, non-submodular functions, optimal experi- mental design, inverse problems, mutual information, uncertainty quantification, Bayesian statistics"
}

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.

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

The estimation of rare event or failure probabilities in high dimensions is of interest in many areas of science and technology. We consider problems where the rare event is expressed in terms of a computationally costly numerical model. Importance sampling with the cross-entropy method offers an efficient way to address such problems provided that a suitable parametric family of biasing densities is employed. Although some existing parametric distribution families are designed to perform efficiently in high dimensions, their applicability within the cross-entropy method is limited to problems with dimension of O(10^2). In this work, rather than directly building sampling densities in high dimensions, we focus on identifying the intrinsic low-dimensional structure of the rare event simulation problem. To this end, we exploit a connection between rare event simulation and Bayesian inverse problems. This allows us to adapt dimension reduction techniques from Bayesian inference to construct new, effectively low-dimensional, biasing distributions within the cross-entropy method. In particular, we employ the approach in [Zahm et al. 2018], as it enables control of the error in the approximation of the optimal biasing distribution. We illustrate our method using two standard high-dimensional reliability benchmark problems and one structural mechanics application involving random fields.

@article{uribe-crossentropy-2020,
  author     = "Uribe, F. and Papaioannou, I. and Marzouk, Y. M. and Straub, D.",
  journal    = "SIAM/ASA Journal on Uncertainty Quantification",
  number     = "2",
  pages      = "818–847",
  title      = "Cross-entropy based importance sampling with failure-informed dimension reduction for rare event simulation",
  volume     = "9",
  year       = "2021",
  doi        = "https://doi.org/10.1137/20M1344585",
  abstract   = "The estimation of rare event or failure probabilities in high dimensions is of interest in many areas of science and technology. We consider problems where the rare event is expressed in terms of a computationally costly numerical model. Importance sampling with the cross-entropy method offers an efficient way to address such problems provided that a suitable parametric family of biasing densities is employed. Although some existing parametric distribution families are designed to perform efficiently in high dimensions, their applicability within the cross-entropy method is limited to problems with dimension of O(10^2). In this work, rather than directly building sampling densities in high dimensions, we focus on identifying the intrinsic low-dimensional structure of the rare event simulation problem. To this end, we exploit a connection between rare event simulation and Bayesian inverse problems. This allows us to adapt dimension reduction techniques from Bayesian inference to construct new, effectively low-dimensional, biasing distributions within the cross-entropy method. In particular, we employ the approach in [Zahm et al. 2018], as it enables control of the error in the approximation of the optimal biasing distribution. We illustrate our method using two standard high-dimensional reliability benchmark problems and one structural mechanics application involving random fields.",
  keywords   = "rare event simulation, reliability analysis, likelihood-informed subspace, importance sampling, cross- entropy method, random fields"
}

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.

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

We propose a framework for solving high-dimensional Bayesian inference problems using \emph{structure-exploiting} low-dimensional transport maps or flows. These maps are confined to a low-dimensional subspace (hence, lazy), and the subspace is identified by minimizing an upper bound on the Kullback–Leibler divergence (hence, structured). Our framework provides a principled way of identifying and exploiting low-dimensional structure in an inference problem. It focuses the expressiveness of a transport map along the directions of most significant discrepancy from the posterior, and can be used to build deep compositions of lazy maps, where low-dimensional projections of the parameters are iteratively transformed to match the posterior. We prove weak convergence of the generated sequence of distributions to the posterior, and we demonstrate the benefits of the framework on challenging inference problems in machine learning and differential equations, using inverse autoregressive flows and polynomial maps as examples of the underlying density estimators.

@article{brennanbigoni-lazy-2020,
  author     = "Brennan, M. and Bigoni, D. and Zahm, O. and Spantini, A. and Marzouk, Y. M.",
  journal    = "Advances in Neural Information Processing Systems (NeurIPS)  oral presentation",
  title      = "Greedy inference with structure-exploiting lazy maps",
  year       = "2020",
  abstract   = "We propose a framework for solving high-dimensional Bayesian inference problems using \emph{structure-exploiting} low-dimensional transport maps or flows. These maps are confined to a low-dimensional subspace (hence, lazy), and the subspace is identified by minimizing an upper bound on the Kullback--Leibler divergence (hence, structured). Our framework provides a principled way of identifying and exploiting low-dimensional structure in an inference problem. It focuses the expressiveness of a transport map along the directions of most significant discrepancy from the posterior, and can be used to build deep compositions of lazy maps, where low-dimensional projections of the parameters are iteratively transformed to match the posterior. We prove weak convergence of the generated sequence of distributions to the posterior, and we demonstrate the benefits of the framework on challenging inference problems in machine learning and differential equations, using inverse autoregressive flows and polynomial maps as examples of the underlying density estimators.",
  keywords   = "Bayesian inference, greedy approximation, transport"
}

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.

@article{bigoni-discretization-2020,
  author     = "Bigoni, D. and Chen, Y. and Garcia-Trillos, N. and Marzouk, Y. M. and Sanz-Alonso, D.",
  journal    = "Inverse Problems",
  number     = "10",
  pages      = "105008",
  title      = "Data-driven forward discretizations for Bayesian inversion",
  volume     = "36",
  year       = "2020",
  doi        = "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."
}

Satellite remote sensing provides a global view to processes on Earth that has unique benefits compared to making measurements on the ground, such as global coverage and enormous data volume. The typical downsides are spatial and temporal gaps and potentially low data quality. Meaningful statistical inference from such data requires overcoming these problems and developing efficient and robust computational tools. We design and implement a computationally efficient multi-scale Gaussian process (GP) software package, satGP, geared towards remote sensing applications. The software is able to handle problems of enormous sizes and to compute marginals and sample from the random field conditioning on at least hundreds of millions of observations. This is achieved by optimizing the computation by, e.g., randomization and splitting the problem into parallel local subproblems which aggressively discard uninformative data.

We describe the mean function of the Gaussian process by approximating marginals of a Markov random field (MRF). Variability around the mean is modeled with a multi-scale covariance kernel, which consists of Matérn, exponential, and periodic components. We also demonstrate how winds can be used to inform covariances locally. The covariance kernel parameters are learned by calculating an approximate marginal maximum likelihood estimate, and the validity of both the multi-scale approach and the method used to learn the kernel parameters is verified in synthetic experiments.

We apply these techniques to a moderate size ozone data set produced by an atmospheric chemistry model and to the very large number of observations retrieved from the Orbiting Carbon Observatory 2 (OCO-2) satellite. The satGP software is released under an open-source license.

@article{susiluoto-satGP-2020,
  author     = "Susiluoto, J. and Spantini, A. and Haario, H. and Härkönen, T. and Marzouk, Y. M.",
  journal    = "Geoscientific Model Development",
  pages      = "3439–3463",
  title      = "Efficient multi-scale Gaussian process regression for massive remote sensing data with satGP v0.1.2.",
  volume     = "13",
  year       = "2020",
  doi        = "10.5194/gmd-13-3439-2020",
  abstract   = "Satellite remote sensing provides a global view to processes on Earth that has unique benefits compared to making measurements on the ground, such as global coverage and enormous data volume. The typical downsides are spatial and temporal gaps and potentially low data quality. Meaningful statistical inference from such data requires overcoming these problems and developing efficient and robust computational tools. We design and implement a computationally efficient multi-scale Gaussian process (GP) software package, satGP, geared towards remote sensing applications. The software is able to handle problems of enormous sizes and to compute marginals and sample from the random field conditioning on at least hundreds of millions of observations. This is achieved by optimizing the computation by, e.g., randomization and splitting the problem into parallel local subproblems which aggressively discard uninformative data.

We describe the mean function of the Gaussian process by approximating marginals of a Markov random field (MRF). Variability around the mean is modeled with a multi-scale covariance kernel, which consists of Matérn, exponential, and periodic components. We also demonstrate how winds can be used to inform covariances locally. The covariance kernel parameters are learned by calculating an approximate marginal maximum likelihood estimate, and the validity of both the multi-scale approach and the method used to learn the kernel parameters is verified in synthetic experiments.

We apply these techniques to a moderate size ozone data set produced by an atmospheric chemistry model and to the very large number of observations retrieved from the Orbiting Carbon Observatory 2 (OCO-2) satellite. The satGP software is released under an open-source license."
}

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.

@article{sciortino-impurity-2020,
  author     = "Sciortino, F. and Howard, N. T. and Marmar, E. S. and Odstrcil, T. and Cao, N. M. and Dux, R. and Hubbard, A. E. and Hughes, J. W. and Irby, J. H. and Marzouk, Y. M. and Milanese, L. M.",
  journal    = "Nuclear Fusion",
  pages      = "126014",
  title      = "Inference of experimental radial impurity transport on Alcator C-Mod: Bayesian parameter estimation and model selection",
  volume     = "60",
  year       = "2020",
  doi        = "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."
}

Multivariate functions encountered in high-dimensional uncertainty quantification problems often vary along a few dominant directions in the input parameter space. We propose a gradient-based method for detecting these directions and using them to construct {ridge approximations} of such functions, in a setting where the functions are vector-valued (e.g., taking values in $\mathbb{R}^n$). The methodology consists of minimizing an upper bound on the approximation error, obtained by {subspace Poincaré inequalities}. We provide a thorough mathematical analysis in the case where the parameter space is equipped with a Gaussian probability measure. The resulting method generalizes the notion of active subspaces associated with scalar-valued functions. A numerical illustration shows that using gradients of the function yields effective dimension reduction. We also show how the choice of norm on the codomain of the function has an impact on the function’s low-dimensional approximation.

@article{zahm-dimred-2020,
  author     = "Zahm, O. and Constantine, P. and Prieur, C. and Marzouk, Y. M.",
  journal    = "SIAM Journal on Scientific Computing",
  number     = "1",
  pages      = "A534-A558",
  title      = "Gradient-based dimension reduction of multivariate vector-valued functions",
  volume     = "42",
  year       = "2020",
  abstract   = "Multivariate functions encountered in high-dimensional uncertainty quantification problems often vary along a few dominant directions in the input parameter space. We propose a gradient-based method for detecting these directions and using them to construct {ridge approximations} of such functions, in a setting where the functions are vector-valued (e.g., taking values in $\mathbb{R}^n$). The methodology consists of minimizing an upper bound on the approximation error, obtained by {subspace Poincaré inequalities}. We provide a thorough mathematical analysis in the case where the parameter space is equipped with a Gaussian probability measure. The resulting method generalizes the notion of active subspaces associated with scalar-valued functions. A numerical illustration shows that using gradients of the function yields effective dimension reduction. We also show how the choice of norm on the codomain of the function has an impact on the function's low-dimensional approximation.",
  keywords   = "high-dimensional function approximation, dimension reduction, active subspace, ridge approximation, Karhunen-Loève decomposition, Poincaré inequality"
}

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.

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

Markov chain Monte Carlo (MCMC) samplers are numerical methods for drawing samples from a given target probability distribution. We discuss one particular MCMC sampler, the MALA-within-Gibbs sampler, from the theoretical and practical perspectives. We first show that the acceptance ratio and step size of this sampler are independent of the overall problem dimension when (i) the target distribution has sparse conditional structure, and (ii) this structure is reflected in the partial updating strategy of MALA-within-Gibbs. If, in addition, the target density is block-wise log-concave, then the sampler’s convergence rate is independent of dimension. From a practical perspective, we expect that MALA-within-Gibbs is useful for solving high-dimensional Bayesian inference problems where the posterior exhibits sparse conditional structure at least approximately. In this context, a partitioning of the state that correctly reflects the sparse conditional structure must be found, and we illustrate this process in two numerical examples. We also discuss trade-offs between the block size used for partial updating and computational requirements that may increase with the number of blocks.

@article{tong-mala-2020,
  author     = "Tong, X. T. and Morzfeld, M. and Marzouk, Y. M.",
  journal    = "SIAM Journal on Scientific Computing",
  number     = "3",
  pages      = "A1765-A1788",
  title      = "MALA-within-Gibbs samplers for high-dimensional distributions with sparse conditional structure",
  volume     = "42",
  year       = "2020",
  abstract   = "Markov chain Monte Carlo (MCMC) samplers are numerical methods for drawing samples from a given target probability distribution. We discuss one particular MCMC sampler, the MALA-within-Gibbs sampler, from the theoretical and practical perspectives. We first show that the acceptance ratio and step size of this sampler are independent of the overall problem dimension when (i) the target distribution has sparse conditional structure, and (ii) this structure is reflected in the partial updating strategy of MALA-within-Gibbs. If, in addition, the target density is block-wise log-concave, then the sampler's convergence rate is independent of dimension. From a practical perspective, we expect that MALA-within-Gibbs is useful for solving high-dimensional Bayesian inference problems where the posterior exhibits sparse conditional structure at least approximately. In this context, a partitioning of the state that correctly reflects the sparse conditional structure must be found, and we illustrate this process in two numerical examples. We also discuss trade-offs between the block size used for partial updating and computational requirements that may increase with the number of blocks.",
  keywords   = "Bayesian computation, high-dimensional distribution, Markov chain Monte Carlo, sparsity, conditional independence, Metropolis-adjusted Langevin"
}

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.

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

Optimization-based samplers such as randomize-then-optimize (RTO) provide an efficient and parallellizable approach to solving large-scale Bayesian inverse problems. These methods solve randomly perturbed optimization problems to draw samples from an approximate posterior distribution. ‘‘Correcting’’ these samples, either by Metropolization or importance sampling, enables characterization of the original posterior distribution. This paper focuses on the scalability of RTO to problems with high- or infinite-dimensional parameters. In particular, we introduce a new subspace strategy to reformulate RTO. For problems with intrinsic low-rank structure, this subspace acceleration makes the computational complexity of RTO scale linearly with the parameter dimension. Furthermore, this subspace perspective suggests a natural extension of RTO to a function space setting. We thus formalize a function space version of RTO and establish sufficient conditions for it to produce a valid Metropolis–Hastings proposal, yielding dimension-independent sampling performance. Numerical examples corroborate the dimension-independence of RTO and demonstrate sampling performance that is also robust to small observational noise.

@article{art_76,
  author     = "Bardsley, J. M. and Cui, T. and Marzouk, Y. M. and Wang, Z.",
  journal    = "SIAM Journal on Scientific Computing",
  number     = "2",
  pages      = "A1317-1347",
  title      = "Scalable optimization-based sampling on function space",
  volume     = "42",
  year       = "2020",
  abstract   = "Optimization-based samplers such as randomize-then-optimize (RTO) provide an efficient and parallellizable approach to solving large-scale Bayesian inverse problems. These methods solve randomly perturbed optimization problems to draw samples from an approximate posterior distribution. ``Correcting'' these samples, either by Metropolization or importance sampling, enables characterization of the original posterior distribution. This paper focuses on the scalability of RTO to problems with high- or infinite-dimensional parameters. In particular, we introduce a new subspace strategy to reformulate RTO. For problems with intrinsic low-rank structure, this subspace acceleration makes the computational complexity of RTO scale linearly with the parameter dimension. Furthermore, this subspace perspective suggests a natural extension of RTO to a function space setting. We thus formalize a function space version of RTO and establish sufficient conditions for it to produce a valid Metropolis--Hastings proposal, yielding dimension-independent sampling performance. Numerical examples corroborate the dimension-independence of RTO and demonstrate sampling performance that is also robust to small observational noise.",
  keywords   = "Markov chain Monte Carlo, Metropolis independence sampling, Bayesian inference, infinite-dimensional inverse problems"
}

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.

@article{art_75,
  author     = "Chilenski, M. and Greenwald, M. and Marzouk, Y. M. and Rice, J. and White, A. E.",
  journal    = "Plasma Physics and Controlled Fusion",
  month      = "11",
  pages      = "125012",
  title      = "On the importance of model selection when inferring impurity transport coefficient profiles",
  volume     = "61",
  year       = "2019",
  doi        = "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."
}

Markov chain Monte Carlo (MCMC) sampling of posterior distributions arising in Bayesian inverse problems is challenging when evaluations of the forward model are computationally expensive. Replacing the forward model with a low-cost, low-fidelity model often significantly reduces computational cost; however, employing a low-fidelity model alone means that the stationary distribution of the MCMC chain is the posterior distribution corresponding to the low-fidelity model, rather than the original posterior distribution corresponding to the high-fidelity model. We propose a multifidelity approach that combines, rather than replaces, the high-fidelity model with a low-fidelity model. First, the low-fidelity model is used to construct a transport map that deterministically couples a reference Gaussian distribution with an approximation of the low-fidelity posterior. Then, the high-fidelity posterior distribution is explored using a non-Gaussian proposal distribution derived from the transport map. This multifidelity "preconditioned" MCMC approach seeks efficient sampling via a proposal that is explicitly tailored to the posterior at hand and that is constructed efficiently with the low-fidelity model. By relying on the low-fidelity model only to construct the proposal distribution, our approach guarantees that the stationary distribution of the MCMC chain is the high-fidelity posterior. In our numerical examples, our multifidelity approach achieves significant speedups compared to single-fidelity MCMC sampling methods.

@article{art_74,
  author     = "Peherstorfer, B. and Marzouk, Y. M.",
  journal    = "Advances in Computational Mathematics",
  month      = "12",
  number     = "5-6",
  pages      = "2321--2348",
  title      = "A transport-based multifidelity preconditioner for Markov chain Monte Carlo",
  volume     = "45",
  year       = "2019",
  doi        = "10.1007/s10444-019-09711-y",
  abstract   = "Markov chain Monte Carlo (MCMC) sampling of posterior distributions arising in Bayesian inverse problems is challenging when evaluations of the forward model are computationally expensive. Replacing the forward model with a low-cost, low-fidelity model often significantly reduces computational cost; however, employing a low-fidelity model alone means that the stationary distribution of the MCMC chain is the posterior distribution corresponding to the low-fidelity model, rather than the original posterior distribution corresponding to the high-fidelity model. We propose a multifidelity approach that combines, rather than replaces, the high-fidelity model with a low-fidelity model. First, the low-fidelity model is used to construct a transport map that deterministically couples a reference Gaussian distribution with an approximation of the low-fidelity posterior. Then, the high-fidelity posterior distribution is explored using a non-Gaussian proposal distribution derived from the transport map. This multifidelity "preconditioned" MCMC approach seeks efficient sampling via a proposal that is explicitly tailored to the posterior at hand and that is constructed efficiently with the low-fidelity model. By relying on the low-fidelity model only to construct the proposal distribution, our approach guarantees that the stationary distribution of the MCMC chain is the high-fidelity posterior. In our numerical examples, our multifidelity approach achieves significant speedups compared to single-fidelity MCMC sampling methods.",
  keywords   = "Bayesian inverse problems, transport maps, multifidelity, model reduction, Markov chain Monte Carlo"
}

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.

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

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

@article{art_67,
  author     = "Morzfeld, M. and Tong, X. T. and Marzouk, Y. M.",
  journal    = "Journal of Computational Physics",
  pages      = "1--28",
  title      = "Localization for MCMC: sampling high-dimensional posterior distributions with local structure",
  volume     = "380",
  year       = "2019",
  doi        = "10.1016/j.jcp.2018.12.008",
  abstract   = "We investigate how ideas from covariance localization in numerical weather prediction can be used in Markov chain Monte Carlo (MCMC) sampling of high-dimensional posterior distributions arising in Bayesian inverse problems. To localize an inverse problem is to enforce an anticipated "local" structure by (i) neglecting small off-diagonal elements of the prior precision and covariance matrices; and (ii) restricting the influence of observations to their neighborhood. For linear problems we can specify the conditions under which posterior moments of the localized problem are close to those of the original problem. We explain physical interpretations of our assumptions about local structure and discuss the notion of high dimensionality in local problems, which is different from the usual notion of high dimensionality in function space MCMC. The Gibbs sampler is a natural choice of MCMC algorithm for localized inverse problems and we demonstrate that its convergence rate is independent of dimension for localized linear problems. Nonlinear problems can also be tackled efficiently by localization and, as a simple illustration of these ideas, we present a localized Metropolis-within-Gibbs sampler. Several linear and nonlinear numerical examples illustrate localization in the context of MCMC samplers for inverse problems.",
  keywords   = "Markov chain Monte Carlo, Bayesian inverse problems, high dimensions, localization, dimension-independent convergence"
}

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.

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

Uncertainty quantification in expensive turbulent combustion simulations usually adopts response surface techniques to accelerate Monte Carlo sampling. However, it is computationally intractable to build response surfaces for high-dimensional kinetic parameters. We employ the active subspaces approach to reduce the dimension of the parameter space, such that building a response surface on the resulting low-dimensional subspace requires many fewer runs of the expensive simulation, rendering the approach suitable for various turbulent combustion models. We demonstrate this approach in simulations of the Cabra H2/N2 jet flame, propagating the uncertainties of 21 kinetic parameters to the liftoff height. We identify a one-dimensional active subspace for the liftoff height using 84 runs of the simulations, from which a response surface with a one-dimensional input is built; the probability distribution of the liftoff height is then characterized by evaluating a large number of samples using the inexpensive response surface. In addition, the active subspace provides a global sensitivity metric for determining the most influential reactions. Comparison with autoignition tests reveals that the sensitivities to the HO2-related reactions in the Cabra flame are promoted by the diffusion processes. The present work demonstrates the capability of active subspaces in quantifying uncertainty in turbulent combustion simulations and provides physical insights into the flame via the active subspace-based sensitivity metric.

@article{art_69,
  author     = "Ji, W. and Ren, Z. and Marzouk, Y. M. and Law, C. K. J. H.",
  journal    = "Proceedings of the Combustion Institute",
  number     = "2",
  pages      = "2175--2182",
  title      = "Quantifying kinetic uncertainty in turbulent combustion simulations using active subspaces",
  volume     = "37",
  year       = "2019",
  doi        = "10.1016/j.proci.2018.06.206",
  abstract   = "Uncertainty quantification in expensive turbulent combustion simulations usually adopts response surface techniques to accelerate Monte Carlo sampling. However, it is computationally intractable to build response surfaces for high-dimensional kinetic parameters. We employ the active subspaces approach to reduce the dimension of the parameter space, such that building a response surface on the resulting low-dimensional subspace requires many fewer runs of the expensive simulation, rendering the approach suitable for various turbulent combustion models. We demonstrate this approach in simulations of the Cabra H2/N2 jet flame, propagating the uncertainties of 21 kinetic parameters to the liftoff height. We identify a one-dimensional active subspace for the liftoff height using 84 runs of the simulations, from which a response surface with a one-dimensional input is built; the probability distribution of the liftoff height is then characterized by evaluating a large number of samples using the inexpensive response surface. In addition, the active subspace provides a global sensitivity metric for determining the most influential reactions. Comparison with autoignition tests reveals that the sensitivities to the HO2-related reactions in the Cabra flame are promoted by the diffusion processes. The present work demonstrates the capability of active subspaces in quantifying uncertainty in turbulent combustion simulations and provides physical insights into the flame via the active subspace-based sensitivity metric.",
  keywords   = "uncertainty quantification, active subspaces, turbulent lifted flames, autoignition"
}

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.

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

We introduce a new framework for efficient sampling from complex probability distributions, using a combination of transport maps and the Metropolis-Hastings rule. The core idea is to use deterministic couplings to transform typical Metropolis proposal mechanisms (e.g., random walks, Langevin methods) into non-Gaussian proposal distributions that can more effectively explore the target density. Our approach adaptively constructs a lower triangular transport map–-an approximation of the Knothe-Rosenblatt rearrangement–-using information from previous Markov chain Monte Carlo (MCMC) states, via the solution of an optimization problem. This optimization problem is convex regardless of the form of the target distribution, and can be solved efficiently without gradient information from the target probability distribution; the target distribution is instead represented via samples. Sequential updates enable efficient and parallelizable adaptation of the map even for large numbers of samples. We show that this approach uses inexact or truncated maps to produce an adaptive MCMC algorithm that is ergodic for the exact target distribution. Numerical demonstrations on a range of parameter inference problems show order-of-magnitude speedups over standard MCMC techniques, measured by the number of effectively independent samples produced per target density evaluation and per unit of wallclock time.

@article{art_66,
  author     = "Parno, M. and Marzouk, Y. M.",
  journal    = "SIAM/ASA Journal on Uncertainty Quantification",
  number     = "2",
  pages      = "645--682",
  title      = "Transport map accelerated Markov chain Monte Carlo",
  volume     = "6",
  year       = "2018",
  doi        = "10.1137/17M1134640",
  abstract   = "We introduce a new framework for efficient sampling from complex probability distributions, using a combination of transport maps and the Metropolis-Hastings rule. The core idea is to use deterministic couplings to transform typical Metropolis proposal mechanisms (e.g., random walks, Langevin methods) into non-Gaussian proposal distributions that can more effectively explore the target density. Our approach adaptively constructs a lower triangular transport map---an approximation of the Knothe-Rosenblatt rearrangement---using information from previous Markov chain Monte Carlo (MCMC) states, via the solution of an optimization problem. This optimization problem is convex regardless of the form of the target distribution, and can be solved efficiently without gradient information from the target probability distribution; the target distribution is instead represented via samples. Sequential updates enable efficient and parallelizable adaptation of the map even for large numbers of samples. We show that this approach uses inexact or truncated maps to produce an adaptive MCMC algorithm that is ergodic for the exact target distribution. Numerical demonstrations on a range of parameter inference problems show order-of-magnitude speedups over standard MCMC techniques, measured by the number of effectively independent samples produced per target density evaluation and per unit of wallclock time.",
  keywords   = "adaptive MCMC, Bayesian inference, measure transformation, optimal transport"
}

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.

@article{art_65,
  author     = "Zhou, Q. and Liu, W. and Li, J. and Marzouk, Y. M.",
  journal    = "Inverse Problems",
  number     = "9",
  pages      = "095001",
  title      = "An approximate empirical Bayesian method for large-scale linear-Gaussian inverse problems",
  volume     = "34",
  year       = "2018",
  doi        = "10.1088/1361-6420/aac287",
  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."
}

This paper presents a methodology for identifying important discipline couplings in multicomponent engineering systems. Coupling among disciplines contributes significantly to the computational cost of analyzing a system, and can become particularly burdensome when coupled analyses are embedded within a design or optimization loop. In many cases, disciplines may be weakly coupled, so that some of the coupling or interaction terms can be neglected without significantly. impacting the accuracy of the system output. Typical practice derives such approximations in an ad hoc manner using expert opinion and domain experience. This work proposes a new approach that formulates an optimization problem to find a model that optimally balances accuracy of the model outputs with the sparsity of the discipline couplings. An adaptive sequential Monte Carlo. sampling-based technique is used to efficiently search the combinatorial model space of different discipline couplings. An algorithm for selecting an optimal model is presented and illustrated in a fire detection satellite model and a turbine engine cycle analysis model.

@article{art_64,
  author     = "Baptista, R. and Marzouk, Y. M. and Willcox, K. and Peherstorfer, B.",
  journal    = "AIAA Journal",
  number     = "6",
  pages      = "2412--2428",
  title      = "Optimal approximations of coupling in multidisciplinary models",
  volume     = "56",
  year       = "2018",
  doi        = "10.2514/1.J056888",
  abstract   = "This paper presents a methodology for identifying important discipline couplings in multicomponent engineering systems. Coupling among disciplines contributes significantly to the computational cost of analyzing a system, and can become particularly burdensome when coupled analyses are embedded within a design or optimization loop. In many cases, disciplines may be weakly coupled, so that some of the coupling or interaction terms can be neglected without significantly. impacting the accuracy of the system output. Typical practice derives such approximations in an ad hoc manner using expert opinion and domain experience.  This work proposes a new approach that formulates an optimization problem to find a model that optimally balances accuracy of the model outputs with the sparsity of the discipline couplings. An adaptive sequential Monte Carlo. sampling-based technique is used to efficiently search the combinatorial model space of different discipline couplings. An algorithm for selecting an optimal model is presented and illustrated in a fire detection satellite model and a turbine engine cycle analysis model."
}

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.

@article{art_63,
  author     = "Mohammadi-Ghazi, R. and Marzouk, Y. M. and Büyüköztürk, O.",
  journal    = "Pattern Recognition",
  pages      = "601--614",
  title      = "Conditional classifiers and boosted conditional Gaussian mixture models for novelty detection",
  volume     = "81",
  year       = "2018",
  doi        = "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.",
  keywords   = "Novelty detection, mixture models, graphical models, conditional dependence, conditional density, additive modeling, boosting, false positive"
}

In this article we develop a new sequential Monte Carlo (SMC) method for multilevel (ML) Monte Carlo estimation. In particular, the method can be used to estimate expectations with respect to a target probability distribution over an infinite-dimensional and non-compact space as given, for example, by a Bayesian inverse problem with Gaussian random field prior. Under suitable assumptions the MLSMC method has the optimal $O(\epsilon^{-2})$ bound on the cost to obtain a mean-square error of $O(\epsilon^{2})$. The algorithm is accelerated by dimension-independent likelihood-informed (DILI) proposals designed for Gaussian priors, leveraging a novel variation which uses empirical sample covariance information in lieu of Hessian information, hence eliminating the requirement for gradient evaluations. The efficiency of the algorithm is illustrated on two examples: inversion of noisy pressure measurements in a PDE model of Darcy flow to recover the posterior distribution of the permeability field, and inversion of noisy measurements of the solution of an SDE to recover the posterior path measure.

@article{art_62,
  author     = "Beskos, A. and Jasra, A. and Law, K. J. H. and Marzouk, Y. M. and Zhou, Y.",
  journal    = "SIAM/ASA Journal on Uncertainty Quantification",
  pages      = "762--786",
  title      = "Multilevel sequential Monte Carlo with dimension-independent likelihood-informed proposals",
  volume     = "6",
  year       = "2018",
  doi        = "10.1137/17M1120993",
  abstract   = "In this article we develop a new sequential Monte Carlo (SMC) method for multilevel (ML) Monte Carlo estimation. In particular, the method can be used to estimate expectations with respect to a target probability distribution over an infinite-dimensional and non-compact space as given, for example, by a Bayesian inverse problem with Gaussian random field prior. Under suitable assumptions the MLSMC method has the optimal $O(\epsilon^{-2})$ bound on the cost to obtain a mean-square error of $O(\epsilon^{2})$. The algorithm is accelerated by dimension-independent likelihood-informed (DILI) proposals designed for Gaussian priors, leveraging a novel variation which uses empirical sample covariance information in lieu of Hessian information, hence eliminating the requirement for gradient evaluations. The efficiency of the algorithm is illustrated on two examples: inversion of noisy pressure measurements in a PDE model of Darcy flow to recover the posterior distribution of the permeability field, and inversion of noisy measurements of the solution of an SDE to recover the posterior path measure.",
  keywords   = "multilevel Monte Carlo, sequential Monte Carlo, Bayesian inverse problem"
}

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.

@article{art_61,
  author     = "Chilenski, M. and Greenwald, M. and Marzouk, Y. M. and Rice, J. and White, A. E.",
  journal    = "Review of Scientific Instruments",
  month      = "1",
  pages      = "013504",
  title      = "Efficient design and verification of diagnostics for impurity transport experiments",
  volume     = "89",
  year       = "2018",
  doi        = "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."
}

Performing Bayesian inference via Markov chain Monte Carlo (MCMC) can be exceedingly expensive when posterior evaluations invoke the evaluation of a. computationally expensive model, such as a system of partial differential equations. In recent work [Conrad et al. JASA 2016, arXiv:1402.1694], we described a framework for constructing and refining local approximations of such models during an MCMC simulation. These posterior–adapted approximations harness regularity of the model to reduce the computational cost of inference while preserving asymptotic exactness of the Markov chain. Here we describe two extensions of that work. First, we prove that samplers running in parallel can collaboratively construct a shared posterior approximation while ensuring ergodicity of each associated chain, providing a novel opportunity for exploiting parallel computation in MCMC. Second, focusing on the Metropolis–adjusted Langevin algorithm, we describe how a proposal distribution can successfully employ gradients and other relevant information extracted from the approximation. We investigate the practical performance of our strategies using two challenging inference problems, the first in subsurface hydrology and the second in glaciology. Using local approximations constructed via parallel chains, we successfully reduce the run time needed to characterize the posterior distributions in these problems from days to hours and from months to days, respectively, dramatically improving the tractability of Bayesian inference.

@article{art_60,
  author     = "Conrad, P. and Davis, A. and Marzouk, Y. M. and Pillai, N. and Smith, A.",
  journal    = "SIAM/ASA Journal on Uncertainty Quantification",
  number     = "1",
  pages      = "339--373",
  title      = "Parallel local approximation MCMC for expensive models",
  volume     = "6",
  year       = "2018",
  doi        = "10.1137/16M1084080",
  abstract   = "Performing Bayesian inference via Markov chain Monte Carlo (MCMC) can be exceedingly expensive when posterior evaluations invoke the evaluation of a. computationally expensive model, such as a system of partial differential equations. In recent work [Conrad et al. JASA 2016, arXiv:1402.1694], we described a framework for constructing and refining local approximations of such models during an MCMC simulation. These posterior--adapted approximations harness regularity of the model to reduce the computational cost of inference while preserving asymptotic exactness of the Markov chain. Here we describe two extensions of that work. First, we prove that samplers running in parallel can collaboratively construct a shared posterior approximation while ensuring ergodicity of each associated chain, providing a novel opportunity for exploiting parallel computation in MCMC. Second, focusing on the Metropolis--adjusted Langevin algorithm, we describe how a proposal distribution can successfully employ gradients and other relevant information extracted from the approximation. We investigate the practical performance of our strategies using two challenging inference problems, the first in subsurface hydrology and the second in glaciology. Using local approximations constructed via parallel chains, we successfully reduce the run time needed to characterize the posterior distributions in these problems from days to hours and from months to days, respectively, dramatically improving the tractability of Bayesian inference.",
  keywords   = "Markov chain Monte Carlo, parallel computing, Metropolis-adjusted Langevin algorithm, Bayesian inference, approximation theory, local regression, surrogate modeling"
}

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.

@article{art_59,
  author     = "Gorodetsky, A. and Karaman, S. and Marzouk, Y. M.",
  journal    = "International Journal of Robotics Research",
  month      = "3",
  number     = "2-3",
  pages      = "340--377",
  title      = "High-dimensional stochastic optimal control using continuous tensor decompositions",
  volume     = "37",
  year       = "2018",
  doi        = "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.",
  keywords   = "stochastic optimal control, motion planning, dynamic programming, tensor decompositions"
}

Characterizing the rheological properties of faults and the evolution of fault friction during seismic slip are fundamental problems in geology and seismology. Recent increases in the frequency of induced earthquakes have intensified the need for robust methods to estimate fault properties. Here, we present a novel approach for aquifer- and fault-property estimation, which combines coupled multiphysics simulation of injection-induced seismicity with adaptive surrogate-based Bayesian inversion. In a synthetic 2D model, we use aquifer pressure, ground displacements, and fault slip measurements during fluid injection to estimate the dynamic fault friction, the critical slip distance, and the aquifer permeability. Our forward model allows us to observe non-monotonic evolutions of shear traction and slip on the fault resulting from the interplay of several physical mechanisms, including injection-induced aquifer expansion, stress transfer along the fault, and slip-induced stress relaxation. This interplay provides the basis for a successful joint inversion of induced seismicity, yielding well-informed Bayesian posterior distributions of dynamic friction and critical slip. We uncover an inverse relationship between dynamic friction and critical slip distance, which is in agreement with the small dynamic friction and large critical slip reported during seismicity on mature faults.

@article{art_58,
  author     = "Jagalur-Mohan, J. and Jha, B. and Wang, Z. and Juanes, R. and Marzouk, Y. M.",
  journal    = "Geophysical Research Letters",
  number     = "3",
  pages      = "1313--1320",
  title      = "Inferring fault frictional and reservoir hydraulic properties from injection-induced seismicity",
  volume     = "45",
  year       = "2018",
  doi        = "10.1002/2017GL075925",
  abstract   = "Characterizing the rheological properties of faults and the evolution of fault friction during seismic slip are fundamental problems in geology and seismology. Recent increases in the frequency of induced earthquakes have intensified the need for robust methods to estimate fault properties. Here, we present a novel approach for aquifer- and fault-property estimation, which combines coupled multiphysics simulation of injection-induced seismicity with adaptive surrogate-based Bayesian inversion. In a synthetic 2D model, we use aquifer pressure, ground displacements, and fault slip measurements during fluid injection to estimate the dynamic fault friction, the critical slip distance, and the aquifer permeability. Our forward model allows us to observe non-monotonic evolutions of shear traction and slip on the fault resulting from the interplay of several physical mechanisms, including injection-induced aquifer expansion, stress transfer along the fault, and slip-induced stress relaxation. This interplay provides the basis for a successful joint inversion of induced seismicity, yielding well-informed Bayesian posterior distributions of dynamic friction and critical slip. We uncover an inverse relationship between dynamic friction and critical slip distance, which is in agreement with the small dynamic friction and large critical slip reported during seismicity on mature faults."
}

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.

@article{art_57,
  author     = "Spantini, A. and Bigoni, D. and Marzouk, Y. M.",
  journal    = "The Journal of Machine Learning Research",
  number     = "66",
  pages      = "1--71",
  title      = "Inference via low-dimensional couplings",
  volume     = "19",
  year       = "2018",
  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.",
  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"
}

Small earthquakes occur due to natural tectonic motions and are induced by oil and gas production processes. In many oil/gas fields and hydrofracking processeshydrofrackings, induced earthquakes result from fluid extraction or injection. The locations and source mechanisms of these earthquakes provide valuable information about the reservoirs. Analysis of induced seismic events has mostly assumed a double-couple source mechanism. However, recent studies have shown a non-negligible percentage of non-double-couple components of source moment tensors in hydraulic fracturing events, assuming a full moment tensor source mechanism. Without uncertainty quantification of the moment tensor solution, it is difficult to determine the reliability of these source models. This study develops a Bayesian method to perform waveform-based full moment tensor inversion and uncertainty quantification for induced seismic events, accounting for both location and velocity model uncertainties. We conduct tests with synthetic events to validate the method, and then apply our newly developed Bayesian inversion approach to real induced seismicity in an oil/gas field in the sultanate of Oman—determining the uncertainties in the source mechanism and in the location of that event.

@article{art_54,
  author     = "Gu, C. and Marzouk, Y. M. and Toksöz, M. N.",
  journal    = "Geophysical Journal International",
  month      = "3",
  number     = "1",
  pages      = "1963--1985",
  title      = "Waveform-based Bayesian full moment tensor inversion and uncertainty determination for induced seismicity in an oil/gas field",
  volume     = "212",
  year       = "2018",
  doi        = "10.1093/gji/ggx517",
  abstract   = "Small earthquakes occur due to natural tectonic motions and are induced by oil and gas production processes. In many oil/gas fields and hydrofracking processeshydrofrackings, induced earthquakes result from fluid extraction or injection. The locations and source mechanisms of these earthquakes provide valuable information about the reservoirs. Analysis of induced seismic events has mostly assumed a double-couple source mechanism. However, recent studies have shown a non-negligible percentage of non-double-couple components of source moment tensors in hydraulic fracturing events, assuming a full moment tensor source mechanism. Without uncertainty quantification of the moment tensor solution, it is difficult to determine the reliability of these source models. This study develops a Bayesian method to perform waveform-based full moment tensor inversion and uncertainty quantification for induced seismic events, accounting for both location and velocity model uncertainties. We conduct tests with synthetic events to validate the method, and then apply our newly developed Bayesian inversion approach to real induced seismicity in an oil/gas field in the sultanate of Oman—determining the uncertainties in the source mechanism and in the location of that event.",
  keywords   = "induced seismicity, Bayesian inversion, full moment tensor, uncertainty quantification"
}

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.

@article{art_53,
  author     = "Marques, A. and Wang, Q. and Marzouk, Y. M.",
  journal    = "AIAA Journal",
  number     = "2",
  pages      = "482--496",
  title      = "Data-driven integral boundary layer modeling for airfoil performance prediction in the laminar regime",
  volume     = "56",
  year       = "2018",
  doi        = "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."
}

Forward propagation of kinetic uncertainty in combustion simulations usually adopts response surface techniques to accelerate Monte Carlo sampling. Yet it is computationally challenging to build response surfaces for high-dimensional input parameters and expensive combustion models. This study uses the active subspace method to identify a low-dimensional subspace of the input space, within which response surfaces can be built. Active subspace methods have previously been developed only for single (scalar) model outputs, however. This paper introduces a new method that can simultaneously approximate the marginal probability density functions of multiple outputs using a single low-dimensional shared subspace. We identify the shared subspace by solving a least-squares system to compute an appropriate combination of single-output active subspaces. Because the identification of the active subspace for each individual output may require a significant number of samples, this process may be computationally intractable for expensive models such as turbulent combustion simulations. Instead, we propose a heuristic approach that learns the relevant subspaces from cheaper combustion models. The performance of the active subspace for a single output, and of the shared subspace for multiple outputs, is first demonstrated with the ignition delay times and laminar flame speeds of hydrogen/air, methane/air, and dimethyl ether (DME)/air mixtures. Then we demonstrate extrapolatory performance of the shared subspace: using a shared subspace trained on the ignition delays at constant volume, we perform forward propagation of kinetic uncertainties through zero-dimensional HCCI simulations—in particular, single-stage ignition of a natural gas/air mixture and two-stage ignition of a DME/air mixture. We show that the shared subspace can accurately reproduce the probability of ignition failure and the probability density of ignition crank angle conditioned on successful ignition, given uncertainty in the kinetics.

@article{art_52,
  author     = "Ji, W. and Wang, J. and Zahm, O. and Marzouk, Y. M. and Yang, B. and Ren, Z. and Law, C. K. J. H.",
  journal    = "Combustion and Flame",
  pages      = "146--157",
  title      = "Shared low-dimensional subspaces for propagating kinetic uncertainty to multiple outputs",
  volume     = "190",
  year       = "2018",
  doi        = "10.1016/j.combustflame.2017.11.021",
  abstract   = "Forward propagation of kinetic uncertainty in combustion simulations usually adopts response surface techniques to accelerate Monte Carlo sampling. Yet it is computationally challenging to build response surfaces for high-dimensional input parameters and expensive combustion models. This study uses the active subspace method to identify a low-dimensional subspace of the input space, within which response surfaces can be built. Active subspace methods have previously been developed only for single (scalar) model outputs, however. This paper introduces a new method that can simultaneously approximate the marginal probability density functions of multiple outputs using a single low-dimensional shared subspace. We identify the shared subspace by solving a least-squares system to compute an appropriate combination of single-output active subspaces. Because the identification of the active subspace for each individual output may require a significant number of samples, this process may be computationally intractable for expensive models such as turbulent combustion simulations. Instead, we propose a heuristic approach that learns the relevant subspaces from cheaper combustion models. The performance of the active subspace for a single output, and of the shared subspace for multiple outputs, is first demonstrated with the ignition delay times and laminar flame speeds of hydrogen/air, methane/air, and dimethyl ether (DME)/air mixtures. Then we demonstrate extrapolatory performance of the shared subspace: using a shared subspace trained on the ignition delays at constant volume, we perform forward propagation of kinetic uncertainties through zero-dimensional HCCI simulations—in particular, single-stage ignition of a natural gas/air mixture and two-stage ignition of a DME/air mixture. We show that the shared subspace can accurately reproduce the probability of ignition failure and the probability density of ignition crank angle conditioned on successful ignition, given uncertainty in the kinetics.",
  keywords   = "Uncertainty propagation, dimension reduction, active subspaces, shared subspaces, multiple outputs"
}

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.

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

It remains an open question to explain the dramatic change in intrinsic rotation induced by slight changes in electron density (White et al. 2013, Phys. Plasmas 20, 056106). One proposed explanation is that momentum transport is sensitive to the second derivatives of the temperature and density profiles (Lee et al. 2015, Plasma Phys. Controlled Fusion 57, 125006), but it is widely considered to be impossible to measure these higher derivatives. In this paper, we show that it is possible to estimate second derivatives of electron density and temperature using a nonparametric regression technique known as Gaussian process regression (GPR). This technique avoids over-constraining the fit by not assuming an explicit functional form for the fitted curve. The uncertainties, obtained rigorously using Markov chain Monte Carlo (MCMC) sampling, are small enough that it is reasonable to explore hypotheses which depend on second derivatives. It is found that the differences in the second derivatives of $n_e$ and $T_e$ between the peaked and hollow rotation cases are rather small, suggesting that changes in the second derivatives are not likely to explain the experimental results.

@article{art_50,
  author     = "Chilenski, M. and Greenwald, M. and Hubbard, A. and Hughes, J. W. and Lee, J. and Marzouk, Y. M. and Rice, J. and White, A. E.",
  journal    = "Nuclear Fusion",
  month      = "9",
  number     = "12",
  pages      = "126013",
  title      = "Experimentally testing the dependence of momentum transport on second derivatives using Gaussian process regression",
  volume     = "57",
  year       = "2017",
  doi        = "10.1088/1741-4326/aa8387",
  abstract   = "It remains an open question to explain the dramatic change in intrinsic rotation induced by slight changes in electron density (White et al. 2013, Phys. Plasmas 20, 056106). One proposed explanation is that momentum transport is sensitive to the second derivatives of the temperature and density profiles (Lee et al. 2015, Plasma Phys. Controlled Fusion 57, 125006), but it is widely considered to be impossible to measure these higher derivatives. In this paper, we show that it is possible to estimate second derivatives of electron density and temperature using a nonparametric regression technique known as Gaussian process regression (GPR). This technique avoids over-constraining the fit by not assuming an explicit functional form for the fitted curve. The uncertainties, obtained rigorously using Markov chain Monte Carlo (MCMC) sampling, are small enough that it is reasonable to explore hypotheses which depend on second derivatives. It is found that the differences in the second derivatives of $n_e$ and $T_e$ between the peaked and hollow rotation cases are rather small, suggesting that changes in the second derivatives are not likely to explain the experimental results."
}

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.

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

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

@article{art_48,
  author     = "Wang, Z. and Bardsley, J. M. and Solonen, A. and Cui, T. and Marzouk, Y. M.",
  journal    = "SIAM Journal on Scientific Computing",
  number     = "5",
  pages      = "S140--S166",
  title      = "Bayesian inverse problems with l1 priors: a randomize-then-optimize approach",
  volume     = "39",
  year       = "2017",
  doi        = "10.1137/16M1080938",
  abstract   = "Prior distributions for Bayesian inference that rely on the $l_1$-norm of the parameters are of considerable interest, in part because they promote parameter fields with less regularity than Gaussian priors (e.g., discontinuities and blockiness). These $l_1$-type priors include the total variation (TV) prior and the Besov space $B^s_{1,1}$ prior, and in general yield non-Gaussian posterior distributions. Sampling from these posteriors is challenging, particularly in the inverse problem setting where the parameter space is high-dimensional and the forward problem may be nonlinear. This paper extends the randomize-then-optimize (RTO) method, an optimization-based sampling algorithm developed for Bayesian inverse problems with Gaussian priors, to inverse problems with $l_1$-type priors. We use a variable transformation to convert an $l_1$-type prior to a standard Gaussian prior, such that the posterior distribution of the transformed parameters is amenable to Metropolized sampling via RTO. We demonstrate this approach on several deconvolution problems and an elliptic PDE inverse problem, using TV or Besov space $B^s_{1,1}$ priors. Our results show that the transformed RTO algorithm characterizes the correct posterior distribution and can be more efficient than other sampling algorithms. The variable transformation can also be extended to other non-Gaussian priors.",
  keywords   = "Inverse problems, Bayesian inference, Monte Carlo methods"
}

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.

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

The accurate approximation of high-dimensional functions is an essential task in uncertainty quantification and many other fields. We propose a new function approximation scheme based on a spectral extension of the tensor-train (TT) decomposition. We first define a functional version of the TT decomposition and analyze its properties. We obtain results on the convergence of the decomposition, revealing links between the regularity of the function, the dimension of the input space, and the TT ranks. We also show that the regularity of the target function is preserved by the univariate functions (i.e., the "cores") comprising the functional TT decomposition. This result motivates an approximation scheme employing polynomial approximations of the cores. For functions with appropriate regularity, the resulting \textit{spectral tensor-train decomposition} combines the favorable dimension-scaling of the TT decomposition with the spectral convergence rate of polynomial approximations, yielding efficient and accurate surrogates for high-dimensional functions. To construct these decompositions, we use the sampling algorithm \texttt{TT-DMRG-cross} to obtain the TT decomposition of tensors resulting from suitable discretizations of the target function. We assess the performance of the method on a range of numerical examples: a modifed set of Genz functions with dimension up to 100, and functions with mixed Fourier modes or with local features. We observe significant improvements in performance over an anisotropic adaptive Smolyak approach. The method is also used to approximate the solution of an elliptic PDE with random input data. The open source software and examples presented in this work are available online.

@article{art_45,
  author     = "Bigoni, D. and Engsig-Karup, A. and Marzouk, Y. M.",
  journal    = "SIAM Journal on Scientific Computing",
  number     = "4",
  pages      = "A2405--A2439",
  title      = "Spectral tensor-train decomposition",
  volume     = "38",
  year       = "2016",
  doi        = "10.1137/15M1036919",
  abstract   = "The accurate approximation of high-dimensional functions is an essential task in uncertainty quantification and many other fields. We propose a new function approximation scheme based on a spectral extension of the tensor-train (TT) decomposition. We first define a functional version of the TT decomposition and analyze its properties. We obtain results on the convergence of the decomposition, revealing links between the regularity of the function, the dimension of the input space, and the TT ranks. We also show that the regularity of the target function is preserved by the univariate functions (i.e., the "cores") comprising the functional TT decomposition. This result motivates an approximation scheme employing polynomial approximations of the cores. For functions with appropriate regularity, the resulting \textit{spectral tensor-train decomposition} combines the favorable dimension-scaling of the TT decomposition with the spectral convergence rate of polynomial approximations, yielding efficient and accurate surrogates for high-dimensional functions. To construct these decompositions, we use the sampling algorithm \texttt{TT-DMRG-cross} to obtain the TT decomposition of tensors resulting from suitable discretizations of the target function. We assess the performance of the method on a range of numerical examples: a modifed set of Genz functions with dimension up to 100, and functions with mixed Fourier modes or with local features. We observe significant improvements in performance over an anisotropic adaptive Smolyak approach. The method is also used to approximate the solution of an elliptic PDE with random input data. The open source software and examples presented in this work are available online.",
  keywords   = "approximation theory; tensor-train decomposition; orthogonal polynomials; uncertainty quantification"
}

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.

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

This paper examines experimental design procedures used to develop surrogates of computational models, exploring the interplay between experimental designs and approximation algorithms. We focus on two widely used approximation approaches, Gaussian process (GP) regression and nonintrusive polynomial approximation. First, we introduce algorithms for minimizing a posterior integrated variance (IVAR) design criterion for GP regression. Our formulation treats design as a continuous optimization problem that can be solved with gradient-based methods on complex input domains without resorting to greedy approximations. We show that minimizing IVAR in this way yields point sets with good interpolation properties and that it enables more accurate GP regression than designs based on entropy minimization or mutual information maximization. Second, using a Mercer kernel/eigenfunction perspective on GP regression, we identify conditions under which GP regression coincides with pseudospectral polynomial approximation. Departures from these conditions can be understood as changes either to the kernel or to the experimental design itself. We then show how IVAR-optimal designs, while sacrificing discrete orthogonality of the kernel eigenfunctions, can yield lower approximation error than orthogonalizing point sets. Finally, we compare the performance of adaptive GP regression and adaptive pseudospectral approximation for several classes of target functions, identifying features that are favorable to the GP + IVAR approach.

@article{art_43,
  author     = "Gorodetsky, A. and Marzouk, Y. M.",
  journal    = "SIAM/ASA Journal on Uncertainty Quantification",
  number     = "1",
  pages      = "796--828",
  title      = "Mercer kernels and integrated variance experimental design: connections between Gaussian process regression and polynomial approximation",
  volume     = "4",
  year       = "2016",
  doi        = "10.1137/15M1017119",
  abstract   = "This paper examines experimental design procedures used to develop surrogates of computational models, exploring the interplay between experimental designs and approximation algorithms. We focus on two widely used approximation approaches, Gaussian process (GP) regression and nonintrusive polynomial approximation. First, we introduce algorithms for minimizing a posterior integrated variance (IVAR) design criterion for GP regression. Our formulation treats design as a continuous optimization problem that can be solved with gradient-based methods on complex input domains without resorting to greedy approximations. We show that minimizing IVAR in this way yields point sets with good interpolation properties and that it enables more accurate GP regression than designs based on entropy minimization or mutual information maximization. Second, using a Mercer kernel/eigenfunction perspective on GP regression, we identify conditions under which GP regression coincides with pseudospectral polynomial approximation. Departures from these conditions can be understood as changes either to the kernel or to the experimental design itself. We then show how IVAR-optimal designs, while sacrificing discrete orthogonality of the kernel eigenfunctions, can yield lower approximation error than orthogonalizing point sets. Finally, we compare the performance of adaptive GP regression and adaptive pseudospectral approximation for several classes of target functions, identifying features that are favorable to the GP + IVAR approach.",
  keywords   = "Gaussian process regression; experimental design; computer experiments; approximation theory; polynomial approximation; kernel interpolation; uncertainty quantification"
}

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.

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

Estimating statistics of model outputs with the Monte Carlo method often requires a large number of model evaluations. This leads to long runtimes if the model is expensive to evaluate. Importance sampling is one approach that can lead to a reduction in the number of model evaluations. Importance sampling uses a biasing distribution to sample the model more efficiently, but generating such a biasing distribution can be difficult and usually also requires model evaluations. A different strategy to speed up Monte Carlo sampling is to replace the computationally expensive high-fidelity model with a computationally cheap surrogate model; however, because the surrogate model outputs are only approximations of the high-fidelity model outputs, the estimate obtained using a surrogate model is in general biased with respect to the estimate obtained using the high-fidelity model. We introduce a multifidelity importance sampling (MFIS) method, which combines evaluations of both the high-fidelity and a surrogate model. It uses a surrogate model to facilitate the construction of the biasing distribution, but relies on a small number of evaluations of the high-fidelity model to derive an unbiased estimate of the statistics of interest. We prove that the MFIS estimate is unbiased even in the absence of accuracy guarantees on the surrogate model itself. The MFIS method can be used with any type of surrogate model, such as projection-based reduced-order models and data-fit models. Furthermore, the MFIS method is applicable to black-box models, i.e., where only inputs and the corresponding outputs of the high-fidelity and the surrogate model are available but not the details of the models themselves. We demonstrate on nonlinear and time-dependent problems that our MFIS method achieves speedups of up to several orders of magnitude compared to Monte Carlo with importance sampling that uses the high-fidelity model only.

@article{art_41,
  author     = "Peherstorfer, B. and Cui, T. and Marzouk, Y. M. and Willcox, K.",
  journal    = "Computer Methods in Applied Mechanics and Engineering",
  pages      = "490--509",
  title      = "Multifidelity importance sampling",
  volume     = "300",
  year       = "2016",
  doi        = "http://dx.doi.org/10.1016/j.cma.2015.12.002",
  abstract   = "Estimating statistics of model outputs with the Monte Carlo method often requires a large number of model evaluations. This leads to long runtimes if the model is expensive to evaluate. Importance sampling is one approach that can lead to a reduction in the number of model evaluations. Importance sampling uses a biasing distribution to sample the model more efficiently, but generating such a biasing distribution can be difficult and usually also requires model evaluations. A different strategy to speed up Monte Carlo sampling is to replace the computationally expensive high-fidelity model with a computationally cheap surrogate model; however, because the surrogate model outputs are only approximations of the high-fidelity model outputs, the estimate obtained using a surrogate model is in general biased with respect to the estimate obtained using the high-fidelity model. We introduce a multifidelity importance sampling (MFIS) method, which combines evaluations of both the high-fidelity and a surrogate model. It uses a surrogate model to facilitate the construction of the biasing distribution, but relies on a small number of evaluations of the high-fidelity model to derive an unbiased estimate of the statistics of interest. We prove that the MFIS estimate is unbiased even in the absence of accuracy guarantees on the surrogate model itself. The MFIS method can be used with any type of surrogate model, such as projection-based reduced-order models and data-fit models. Furthermore, the MFIS method is applicable to black-box models, i.e., where only inputs and the corresponding outputs of the high-fidelity and the surrogate model are available but not the details of the models themselves. We demonstrate on nonlinear and time-dependent problems that our MFIS method achieves speedups of up to several orders of magnitude compared to Monte Carlo with importance sampling that uses the high-fidelity model only.",
  keywords   = "Monte Carlo; importance sampling; surrogate modeling; multifidelity methods"
}

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.

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

We construct a new framework for accelerating Markov chain Monte Carlo in posterior sampling problems where standard methods are limited by the computational cost of the likelihood, or of numerical models embedded therein. Our approach introduces local approximations of these models into the Metropolis-Hastings kernel, borrowing ideas from deterministic approximation theory, optimization, and experimental design. Previous efforts at integrating approximate models into inference typically sacrifice either the sampler’s exactness or efficiency; our work seeks to address these limitations by exploiting useful convergence characteristics of local approximations. We prove the ergodicity of our approximate Markov chain, showing that it samples asymptotically from the \emph{exact} posterior distribution of interest. We describe variations of the algorithm that employ either local polynomial approximations or local Gaussian process regressors. Our theoretical results reinforce the key observation underlying this paper: when the likelihood has some \emph{local} regularity, the number of model evaluations per MCMC step can be greatly reduced without biasing the Monte Carlo average. Numerical experiments demonstrate multiple order-of-magnitude reductions in the number of forward model evaluations used in representative ODE and PDE inference problems, with both synthetic and real data.

@article{art_36,
  author     = "Conrad, P. and Marzouk, Y. M. and Pillai, N. and Smith, A.",
  journal    = "Journal of the American Statistical Association",
  number     = "516",
  pages      = "1591--1607",
  title      = "Accelerating asymptotically exact MCMC for computationally intensive models via local approximations",
  volume     = "111",
  year       = "2016",
  doi        = "10.1080/01621459.2015.1096787",
  abstract   = "We construct a new framework for accelerating Markov chain Monte Carlo in posterior sampling problems where standard methods are limited by the computational cost of the likelihood, or of numerical models embedded therein. Our approach introduces local approximations of these models into the Metropolis-Hastings kernel, borrowing ideas from deterministic approximation theory, optimization, and experimental design. Previous efforts at integrating approximate models into inference typically sacrifice either the sampler's exactness or efficiency; our work seeks to address these limitations by exploiting useful convergence characteristics of local approximations. We prove the ergodicity of our approximate Markov chain, showing that it samples asymptotically from the \emph{exact} posterior distribution of interest. We describe variations of the algorithm that employ either local polynomial approximations or local Gaussian process regressors. Our theoretical results reinforce the key observation underlying this paper: when the likelihood has some \emph{local} regularity, the number of model evaluations per MCMC step can be greatly reduced without biasing the Monte Carlo average. Numerical experiments demonstrate multiple order-of-magnitude reductions in the number of forward model evaluations used in representative ODE and PDE inference problems, with both synthetic and real data.",
  keywords   = "Markov chain Monte Carlo, experimental design, approximation theory, local approximation, computer experiments, emulators"
}

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.

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

Process variations can significantly degrade device performance and chip yield in silicon photonics. In order to reduce the design and production costs, it is highly desirable to predict the statistical behavior of a device before the final fabrication. Monte Carlo is the mainstream computational technique used to estimate the uncertainties caused by process variations. However, it is very often too expensive due to its slow convergence rate. Recently, stochastic spectral methods based on polynomial chaos expansions have emerged as a promising alternative, and they have shown significant speedup over Monte Carlo in many engineering problems. The existing literature mostly assumes that the random parameters are mutually independent. However, in practical applications such assumption may not be necessarily accurate. In this paper, we develop an efficient numerical technique based on stochastic collocation to simulate silicon photonics with correlated and non-Gaussian random parameters. The effectiveness of our proposed technique is demonstrated by the simulation results of a silicon-on-insulator based directional coupler example. Since the mathematic formulation in this paper is very generic, our proposed algorithm can be applied to a large class of photonic design cases as well as to many other engineering problems.

@article{art_34,
  author     = "Weng, T. and Zhang, Z. and Su, Z. and Marzouk, Y. M. and Melloni, A. and Daniel, L.",
  journal    = "Optics Express",
  number     = "4",
  pages      = "4242--4254",
  title      = "Uncertainty quantification of silicon photonic devices with correlated and non-Gaussian random parameters",
  volume     = "23",
  year       = "2015",
  doi        = "10.1364/OE.23.004242",
  abstract   = "Process variations can significantly degrade device performance and chip yield in silicon photonics. In order to reduce the design and production costs, it is highly desirable to predict the statistical behavior of a device before the final fabrication. Monte Carlo is the mainstream computational technique used to estimate the uncertainties caused by process variations. However, it is very often too expensive due to its slow convergence rate. Recently, stochastic spectral methods based on polynomial chaos expansions have emerged as a promising alternative, and they have shown significant speedup over Monte Carlo in many engineering problems. The existing literature mostly assumes that the random parameters are mutually independent. However, in practical applications such assumption may not be necessarily accurate. In this paper, we develop an efficient numerical technique based on stochastic collocation to simulate silicon photonics with correlated and non-Gaussian random parameters. The effectiveness of our proposed technique is demonstrated by the simulation results of a silicon-on-insulator based directional coupler example. Since the mathematic formulation in this paper is very generic, our proposed algorithm can be applied to a large class of photonic design cases as well as to many other engineering problems."
}

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.

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

The need to fit smooth temperature and density profiles to discrete observations is ubiquitous in plasma physics, but the prevailing techniques for this have many shortcomings that cast doubt on the statistical validity of the results. This issue is amplified in the context of validation of gyrokinetic transport models (Holland et al. 2009, Phys. Plasmas 16, 052301), where the strong sensitivity of the code outputs to input gradients means that inadequacies in the profile fitting technique can easily lead to an incorrect assessment of the degree of agreement with experimental measurements. In order to rectify the shortcomings of standard approaches to profile fitting, we have applied Gaussian process regression (GPR), a powerful nonparametric regression technique, to analyze an Alcator C-Mod L-mode discharge used for past gyrokinetic validation work (Howard et al. 2012, Nucl. Fusion 52, 063002). We show that the GPR techniques can reproduce the previous results while delivering more statistically rigorous fits and uncertainty estimates for both the value and the gradient of plasma profiles with an improved level of automation. We also discuss how the use of GPR can allow for dramatic increases in the rate of convergence of uncertainty propagation for any code that takes experimental profiles as inputs. The new GPR techniques for profile fitting and uncertainty propagation are quite useful and general, and we describe the steps to implementation in detail in this paper. These techniques have the potential to substantially improve the quality of uncertainty estimates on profile fits and the rate of convergence of uncertainty propagation, making them of great interest for wider use in fusion experiments and modeling efforts.

@article{art_32,
  author     = "Chilenski, M. A. and Greenwald, M. and Marzouk, Y. M. and Howard, N. T. and White, A. E. and Rice, J. and Walk, J. R.",
  journal    = "Nuclear Fusion",
  number     = "2",
  pages      = "023012",
  title      = "Improved profile fitting and quantification of uncertainty in experimental measurements of impurity transport coefficients using Gaussian process regression",
  volume     = "55",
  year       = "2015",
  doi        = "10.1088/0029-5515/55/2/023012",
  abstract   = "The need to fit smooth temperature and density profiles to discrete observations is ubiquitous in plasma physics, but the prevailing techniques for this have many shortcomings that cast doubt on the statistical validity of the results. This issue is amplified in the context of validation of gyrokinetic transport models (Holland et al. 2009, Phys. Plasmas 16, 052301), where the strong sensitivity of the code outputs to input gradients means that inadequacies in the profile fitting technique can easily lead to an incorrect assessment of the degree of agreement with experimental measurements. In order to rectify the shortcomings of standard approaches to profile fitting, we have applied Gaussian process regression (GPR), a powerful nonparametric regression technique, to analyze an Alcator C-Mod L-mode discharge used for past gyrokinetic validation work (Howard et al. 2012, Nucl. Fusion 52, 063002). We show that the GPR techniques can reproduce the previous results while delivering more statistically rigorous fits and uncertainty estimates for both the value and the gradient of plasma profiles with an improved level of automation. We also discuss how the use of GPR can allow for dramatic increases in the rate of convergence of uncertainty propagation for any code that takes experimental profiles as inputs. The new GPR techniques for profile fitting and uncertainty propagation are quite useful and general, and we describe the steps to implementation in detail in this paper. These techniques have the potential to substantially improve the quality of uncertainty estimates on profile fits and the rate of convergence of uncertainty propagation, making them of great interest for wider use in fusion experiments and modeling efforts.",
  keywords   = "plasma physics, Gaussian process regression, uncertainty propagation, model validation"
}

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.

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

Bayesian inference provides a natural framework for combining experimental data with prior knowledge to develop chemical kinetic models and quantify the associated uncertainties, not only in parameter values but also in model structure. Most existing applications of Bayesian model selection methods to chemical kinetics have been limited to comparisons among a small set of models, however. The significant computational cost of evaluating posterior model probabilities renders traditional Bayesian methods infeasible when the model space becomes large. We present a new framework for tractable Bayesian model inference and uncertainty quantification using a large number of systematically generated model hypotheses. The approach involves imposing point-mass mixture priors over rate constants and exploring the resulting posterior distribution using an adaptive Markov chain Monte Carlo method. The posterior samples are used to identify plausible models, to quantify rate constant uncertainties, and to extract key diagnostic information about model structure–-such as the reactions and operating pathways most strongly supported by the data. We provide numerical demonstrations of the proposed framework by inferring kinetic models for catalytic steam and dry reforming of methane using available experimental data.

@article{art_30,
  author     = "Galagali, N. and Marzouk, Y. M.",
  journal    = "Chemical Engineering Science",
  month      = "2",
  pages      = "170--190",
  title      = "Bayesian inference of chemical kinetic models from proposed reactions",
  volume     = "123",
  year       = "2015",
  doi        = "doi:10.1016/j.ces.2014.10.030",
  abstract   = "Bayesian inference provides a natural framework for combining experimental data with prior knowledge to develop chemical kinetic models and quantify the associated uncertainties, not only in parameter values but also in model structure. Most existing applications of Bayesian model selection methods to chemical kinetics have been limited to comparisons among a small set of models, however. The significant computational cost of evaluating posterior model probabilities renders traditional Bayesian methods infeasible when the model space becomes large. We present a new framework for tractable Bayesian model inference and uncertainty quantification using a large number of systematically generated model hypotheses. The approach involves imposing point-mass mixture priors over rate constants and exploring the resulting posterior distribution using an adaptive Markov chain Monte Carlo method. The posterior samples are used to identify plausible models, to quantify rate constant uncertainties, and to extract key diagnostic information about model structure---such as the reactions and operating pathways most strongly supported by the data. We provide numerical demonstrations of the proposed framework by inferring kinetic models for catalytic steam and dry reforming of methane using available experimental data.",
  keywords   = "Bayesian inference, chemical kinetics, model selection, Markov chain Monte Carlo, adaptive MCMC, online expectation maximization"
}

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.

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

The intrinsic dimensionality of an inverse problem is affected by prior information, the accuracy and number of observations, and the smoothing properties of the forward operator. From a Bayesian perspective, changes from the prior to the posterior may, in many problems, be confined to a relatively low-dimensional subspace of the parameter space. We present a dimension reduction approach that defines and identifies such a subspace, called the "likelihood-informed subspace" (LIS), by characterizing the relative influences of the prior and the likelihood over the support of the posterior distribution. This identification enables new and more efficient computational methods for Bayesian inference with nonlinear forward models and Gaussian priors. In particular, we approximate the posterior distribution as the product of a lower-dimensional posterior defined on the LIS and the prior distribution marginalized onto the complementary subspace. Markov chain Monte Carlo sampling can then proceed in lower dimensions, with significant gains in computational efficiency. We also introduce a Rao-Blackwellization strategy that de-randomizes Monte Carlo estimates of posterior expectations for additional variance reduction. We demonstrate the efficiency of our methods using two numerical examples: inference of permeability in a groundwater system governed by an elliptic PDE, and an atmospheric remote sensing problem based on Global Ozone Monitoring System (GOMOS) observations.

@article{art_28,
  author     = "Cui, T. and Martin, J. and Marzouk, Y. M. and Solonen, A. and Spantini, A.",
  journal    = "Inverse Problems",
  pages      = "114015",
  title      = "Likelihood-informed dimension reduction for nonlinear inverse problems",
  volume     = "29",
  year       = "2014",
  doi        = "doi:10.1088/0266-5611/30/11/114015",
  abstract   = "The intrinsic dimensionality of an inverse problem is affected by prior information, the accuracy and number of observations, and the smoothing properties of the forward operator. From a Bayesian perspective, changes from the prior to the posterior may, in many problems, be confined to a relatively low-dimensional subspace of the parameter space. We present a dimension reduction approach that defines and identifies such a subspace, called the "likelihood-informed subspace" (LIS), by characterizing the relative influences of the prior and the likelihood over the support of the posterior distribution. This identification enables new and more efficient computational methods for Bayesian inference with nonlinear forward models and Gaussian priors. In particular, we approximate the posterior distribution as the product of a lower-dimensional posterior defined on the LIS and the prior distribution marginalized onto the complementary subspace. Markov chain Monte Carlo sampling can then proceed in lower dimensions, with significant gains in computational efficiency. We also introduce a Rao-Blackwellization strategy that de-randomizes Monte Carlo estimates of posterior expectations for additional variance reduction. We demonstrate the efficiency of our methods using two numerical examples: inference of permeability in a groundwater system governed by an elliptic PDE, and an atmospheric remote sensing problem based on Global Ozone Monitoring System (GOMOS) observations.",
  keywords   = "Inverse problem, Bayesian inference, dimension reduction, low-rank approximation, Markov chain Monte Carlo, variance reduction"
}

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.

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

The Bayesian approach to inverse problems typically relies on posterior sampling approaches, such as Markov chain Monte Carlo, for which the generation of each sample requires one or more evaluations of the parameter-to-observable map or forward model. When these evaluations are computationally intensive, approximations of the forward model are essential to accelerating sample-based inference. Yet the construction of globally accurate approximations for nonlinear forward models can be computationally prohibitive and in fact unnecessary, as the posterior distribution typically concentrates on a small fraction of the support of the prior distribution. We present a new approach that uses stochastic optimization to construct polynomial approximations over a sequence of measures adaptively determined from the data, eventually concentrating on the posterior distribution. The approach yields substantial gains in efficiency and accuracy over prior-based surrogates, as demonstrated via application to inverse problems in partial differential equations.

@article{art_26,
  author     = "Li, J. and Marzouk, Y. M.",
  journal    = "SIAM Journal on Scientific Computing",
  number     = "3",
  pages      = "A1163--A1186",
  title      = "Adaptive construction of surrogates for the Bayesian solution of inverse problems",
  volume     = "36",
  year       = "2014",
  abstract   = "The Bayesian approach to inverse problems typically relies on posterior sampling approaches, such as Markov chain Monte Carlo, for which the generation of each sample requires one or more evaluations of the parameter-to-observable map or forward model. When these evaluations are computationally intensive, approximations of the forward model are essential to accelerating sample-based inference. Yet the construction of globally accurate approximations for nonlinear forward models can be computationally prohibitive and in fact unnecessary, as the posterior distribution typically concentrates on a small fraction of the support of the prior distribution. We present a new approach that uses stochastic optimization to construct polynomial approximations over a sequence of measures adaptively determined from the data, eventually concentrating on the posterior distribution. The approach yields substantial gains in efficiency and accuracy over prior-based surrogates, as demonstrated via application to inverse problems in partial differential equations.",
  keywords   = "Bayesian inference, cross-entropy method, importance sampling, inverse problem, Kullback-Leibler divergence, Markov chain Monte Carlo, polynomial chaos"
}

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.

@article{art_25,
  author     = "Huan, X. and Marzouk, Y. M.",
  journal    = "International Journal for Uncertainty Quantification",
  number     = "6",
  pages      = "479--510",
  title      = "Gradient-based stochastic optimization methods in Bayesian experimental design",
  volume     = "4",
  year       = "2014",
  doi        = "10.1615/Int.J.UncertaintyQuantification.2014006730",
  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."
}

This work explores the implementation of an adaptive strategy to design sparse ensembles of oceanic simulations suitable for constructing polynomial chaos surrogates. We use a recently developed pseudo-spectral algorithm that is based on a direct application of the Smolyak sparse grid formula and that allows the use of arbitrary admissible sparse grids. The adaptive algorithm is tested using an existing simulation database of the oceanic response to Hurricane Ivan in the Gulf of Mex- ico. The a priori tests demonstrate that sparse and adaptive pseudo-spectral constructions lead to substantial savings over isotropic sparse sampling in the present setting.

@article{art_24,
  author     = "Winokur, J. and Conrad, P. and Sraj, I. and Knio, O. and Srinivasan, A. and Thacker, W. C. and Marzouk, Y. M. and Iskandarani, M.",
  journal    = "Computational Geosciences",
  number     = "6",
  pages      = "899-911",
  title      = "A priori testing of sparse adaptive polynomial chaos expansions using an ocean general circulation model database",
  volume     = "17",
  year       = "2013",
  doi        = "10.1007/s10596-013-9361-3",
  abstract   = "This work explores the implementation of an adaptive strategy to design sparse ensembles of oceanic simulations suitable for constructing polynomial chaos surrogates. We use a recently developed pseudo-spectral algorithm that is based on a direct application of the Smolyak sparse grid formula and that allows the use of arbitrary admissible sparse grids. The adaptive algorithm is tested using an existing simulation database of the oceanic response to Hurricane Ivan in the Gulf of Mex- ico. The a priori tests demonstrate that sparse and adaptive pseudo-spectral constructions lead to substantial savings over isotropic sparse sampling in the present setting."
}

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.

@article{art_23,
  author     = "Conrad, P. and Marzouk, Y. M.",
  journal    = "SIAM Journal on Scientific Computing",
  number     = "6",
  pages      = "A2643--A2670",
  title      = "Adaptive Smolyak pseudospectral approximations",
  volume     = "35",
  year       = "2013",
  doi        = "10.1137/120890715",
  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.",
  keywords   = "Smolyak algorithms, sparse grids, orthogonal polynomials, pseudospectral approximation, approximation theory, uncertainty quantification"
}

The full application of Bayesian inference to inverse problems requires exploration of a posterior distribution that typically does not possess a standard form. In this context, Markov chain Monte Carlo (MCMC) methods are often used. These methods require many evaluations of a computationally intensive forward model to produce the equivalent of one independent sample from the posterior. We consider applications in which approximate forward models at multiple resolution levels are available, each endowed with a probabilistic error estimate. These situations occur, for example, when the forward model involves Monte Carlo integration. We present a novel MCMC method called $MC^3$ that uses low-resolution forward models to approximate draws from a posterior distribution built with the high-resolution forward model. The acceptance ratio is estimated with some statistical error; then a confidence interval for the true acceptance ratio is found, and acceptance is performed correctly with some confidence. The high-resolution models are rarely run and a significant speed up is achieved. Our multiple-resolution forward models themselves are built around a new importance sampling scheme that allows Monte Carlo forward models to be used efficiently in inverse problems. The method is used to solve an inverse transport problem that finds applications in atmospheric remote sensing. We present a path-recycling methodology to efficiently vary parameters in the transport equation. The forward transport equation is solved by a Monte Carlo method that is amenable to the use of $MC^3$ to solve the inverse transport problem using a Bayesian formalism.

@article{art_22,
  author     = "Bal, G. and Langmore, I. and Marzouk, Y. M.",
  journal    = "Inverse Problems and Imaging",
  number     = "1",
  pages      = "81--105",
  title      = "Bayesian inverse problems with Monte Carlo forward models",
  volume     = "7",
  year       = "2013",
  doi        = "10.3934/ipi.2013.7.81",
  abstract   = "The full application of Bayesian inference to inverse problems requires exploration of a posterior distribution that typically does not possess a standard form. In this context, Markov chain Monte Carlo (MCMC) methods are often used. These methods require many evaluations of a computationally intensive forward model to produce the equivalent of one independent sample from the posterior. We consider applications in which approximate forward models at multiple resolution levels are available, each endowed with a probabilistic error estimate. These situations occur, for example, when the forward model involves Monte Carlo integration. We present a novel MCMC method called $MC^3$ that uses low-resolution forward models to approximate draws from a posterior distribution built with the high-resolution forward model. The acceptance ratio is estimated with some statistical error; then a confidence interval for the true acceptance ratio is found, and acceptance is performed correctly with some confidence. The high-resolution models are rarely run and a significant speed up is achieved. Our multiple-resolution forward models themselves are built around a new importance sampling scheme that allows Monte Carlo forward models to be used efficiently in inverse problems. The method is used to solve an inverse transport problem that finds applications in atmospheric remote sensing. We present a path-recycling methodology to efficiently vary parameters in the transport equation. The forward transport equation is solved by a Monte Carlo method that is amenable to the use of $MC^3$ to solve the inverse transport problem using a Bayesian formalism.",
  keywords   = "linear transport; perturbation Monte Carlo; Bayesian; importance sampling; inverse problems; Markov chain Monte Carlo"
}

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.

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

We present a new approach to Bayesian inference that entirely avoids Markov chain simulation, by constructing a map that pushes forward the prior measure to the posterior measure. Existence and uniqueness of a suitable measure-preserving map is established by formulating the problem in the context of optimal transport theory. We discuss various means of explicitly parameterizing the map and computing it efficiently through solution of an optimization problem, exploiting gradient information from the forward model when possible. The resulting algorithm overcomes many of the computational bottlenecks associated with Markov chain Monte Carlo. Advantages of a map-based representation of the posterior include analytical expressions for posterior moments and the ability to generate arbitrary numbers of independent posterior samples without additional likelihood evaluations or forward solves. The optimization approach also provides clear convergence criteria for posterior approximation and facilitates model selection through automatic evaluation of the marginal likelihood. We demonstrate the accuracy and efficiency of the approach on nonlinear inverse problems of varying dimension, involving the inference of parameters appearing in ordinary and partial differential equations.

@article{art_20,
  author     = "Moselhy, T. A. and Marzouk, Y. M.",
  journal    = "Journal of Computational Physics",
  number     = "23",
  pages      = "7815--7850",
  title      = "Bayesian inference with optimal maps",
  volume     = "231",
  year       = "2012",
  doi        = "10.1016/j.jcp.2012.07.022",
  abstract   = "We present a new approach to Bayesian inference that entirely avoids Markov chain simulation, by constructing a map that pushes forward the prior measure to the posterior measure. Existence and uniqueness of a suitable measure-preserving map is established by formulating the problem in the context of optimal transport theory. We discuss various means of explicitly parameterizing the map and computing it efficiently through solution of an optimization problem, exploiting gradient information from the forward model when possible. The resulting algorithm overcomes many of the computational bottlenecks associated with Markov chain Monte Carlo. Advantages of a map-based representation of the posterior include analytical expressions for posterior moments and the ability to generate arbitrary numbers of independent posterior samples without additional likelihood evaluations or forward solves. The optimization approach also provides clear convergence criteria for posterior approximation and facilitates model selection through automatic evaluation of the marginal likelihood. We demonstrate the accuracy and efficiency of the approach on nonlinear inverse problems of varying dimension, involving the inference of parameters appearing in ordinary and partial differential equations.",
  keywords   = "Bayesian inference; Optimal transport; Measure-preserving maps; Inverse problems; Polynomial chaos; Numerical optimization"
}

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

Data assimilation is an essential tool for predicting the behavior of real physical systems given approximate simulation models and limited observations. For many complex systems, there may exist several models, each with different properties and predictive capabilities. It is desirable to incorporate multiple models into the assimilation procedure in order to obtain a more accurate prediction of the physics than any model alone can provide. In this paper, we propose a framework for conducting sequential data assimilation with multiple models and sources of data. The assimilated solution is a linear combination of all model predictions and data. One notable feature is that the combination takes the most general form with matrix weights. By doing so the method can readily utilize different weights in different sections of the solution state vectors, allow the models and data to have different dimensions, and deal with the case of a singular state covariance. We prove that the proposed assimilation method, termed direct assimilation, minimizes a variational functional, a generalized version of the one used in the classical Kalman filter. We also propose an efficient iterative assimilation method that assimilates two models at a time until all models and data are assimilated. The mathematical equivalence of the iterative method and the direct method is established. Numerical examples are presented to demonstrate the effectiveness of the new method.

@article{art_18,
  author     = "Narayan, A. and Marzouk, Y. M. and Xiu, D.",
  journal    = "Journal of Computational Physics",
  number     = "19",
  pages      = "6401--6418",
  title      = "Sequential data assimilation with multiple models",
  volume     = "231",
  year       = "2012",
  doi        = "10.1016/j.jcp.2012.06.002",
  abstract   = "Data assimilation is an essential tool for predicting the behavior of real physical systems given approximate simulation models and limited observations. For many complex systems, there may exist several models, each with different properties and predictive capabilities. It is desirable to incorporate multiple models into the assimilation procedure in order to obtain a more accurate prediction of the physics than any model alone can provide. In this paper, we propose a framework for conducting sequential data assimilation with multiple models and sources of data. The assimilated solution is a linear combination of all model predictions and data. One notable feature is that the combination takes the most general form with matrix weights. By doing so the method can readily utilize different weights in different sections of the solution state vectors, allow the models and data to have different dimensions, and deal with the case of a singular state covariance. We prove that the proposed assimilation method, termed direct assimilation, minimizes a variational functional, a generalized version of the one used in the classical Kalman filter. We also propose an efficient iterative assimilation method that assimilates two models at a time until all models and data are assimilated. The mathematical equivalence of the iterative method and the direct method is established. Numerical examples are presented to demonstrate the effectiveness of the new method.",
  keywords   = "Uncertainty quantification; Data assimilation; Kalman filter; Model averaging"
}

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.

@article{art_17,
  author     = "Berry, R. D. and Najm, H. N. and Debusschere, B. J. and Marzouk, Y. M. and Adalsteinsson, H.",
  journal    = "Journal of Computational Physics",
  number     = "5",
  pages      = "2180-2198",
  title      = "Data-free inference of the joint distribution of uncertain model parameters",
  volume     = "231",
  year       = "2012",
  doi        = "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.",
  keywords   = "Uncertainty quantification; Bayesian statistics; Missing information"
}

This paper presents a novel tabulation strategy for the adaptive numerical integration of chemical kinetics using the computational singular perturbation (CSP) method. The strategy stores and reuses CSP quantities required to filter out fast dissipative processes, resulting in a non-stiff chemical source term. In particular, non-parametric regression on low-dimensional slow invariant manifolds (SIMs) in the chemical state space is used to approximate the CSP vectors spanning the fast chemical subspace and the associated fast chemical time-scales. The relevant manifold and its dimension varies depending on the local number of exhausted modes at every location in the chemical state space. Multiple manifolds are therefore tabulated, corresponding to different numbers of exhausted modes (dimensions) and associated radical species. Non-parametric representations are inherently adaptive, and rely on efficient approximate-nearest-neighbor queries. As the CSP information is only a function of the non-radical species in the system and has relatively small gradients in the chemical state space, tabulation occurs in a lower-dimensional state space and at a relatively coarse level, thereby improving scalability to larger chemical mechanisms. The approach is demonstrated on the simulation of homogeneous constant pressure H2–air and CH4–air ignition, over a range of initial conditions. For CH4–air, results are shown that outperform direct implicit integration of the stiff chemical kinetics while maintaining good accuracy.

@article{art_16,
  author     = "Debusschere, B. J. and Marzouk, Y. M. and Najm, H. N. and Rhoads, B. and Goussis, D. A. and Valorani, M.",
  journal    = "Combustion Theory and Modeling",
  number     = "1",
  pages      = "173-198",
  title      = "Computational singular perturbation with non-parametric tabulation of slow manifolds for time integration of stiff chemical kinetics",
  volume     = "16",
  year       = "2011",
  doi        = "10.1080/13647830.2011.596575",
  abstract   = "This paper presents a novel tabulation strategy for the adaptive numerical integration of chemical kinetics using the computational singular perturbation (CSP) method. The strategy stores and reuses CSP quantities required to filter out fast dissipative processes, resulting in a non-stiff chemical source term. In particular, non-parametric regression on low-dimensional slow invariant manifolds (SIMs) in the chemical state space is used to approximate the CSP vectors spanning the fast chemical subspace and the associated fast chemical time-scales. The relevant manifold and its dimension varies depending on the local number of exhausted modes at every location in the chemical state space. Multiple manifolds are therefore tabulated, corresponding to different numbers of exhausted modes (dimensions) and associated radical species. Non-parametric representations are inherently adaptive, and rely on efficient approximate-nearest-neighbor queries. As the CSP information is only a function of the non-radical species in the system and has relatively small gradients in the chemical state space, tabulation occurs in a lower-dimensional state space and at a relatively coarse level, thereby improving scalability to larger chemical mechanisms. The approach is demonstrated on the simulation of homogeneous constant pressure H2--air and CH4--air ignition, over a range of initial conditions. For CH4--air, results are shown that outperform direct implicit integration of the stiff chemical kinetics while maintaining good accuracy.",
  keywords   = "chemical kinetics, computational singular perturbation, slow manifold, non-parametric regression, nearest neighbors, kd-trees"
}

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.

@article{art_15,
  author     = "McKenna, S. and Ray, J. and Marzouk, Y. M. and van Bloemen Waanders, B.",
  journal    = "Advances in Water Resources",
  month      = "5",
  number     = "5",
  pages      = "617-626",
  title      = "Truncated multi-Gaussian fields and effective conductance of binary media",
  volume     = "34",
  year       = "2011",
  doi        = "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.",
  keywords   = "Upscaling; Binary media; Effective conductivity"
}

Using high-resolution 3-D vortex simulations, this study seeks a mechanistic understanding of vorticity dynamics in transverse jets at a finite Reynolds number. A full no-slip boundary condition, rigorously formulated in terms of vorticity generation along the channel wall, captures unsteady interactions between the wall boundary layer and the jet – in particular, the separation of the wall boundary layer and its transport into the interior. For comparison, we also implement a reduced boundary condition that suppresses the separation of the wall boundary layer away from the jet nozzle. By contrasting results obtained with these two boundary conditions, we characterize near-field vortical structures formed as the wall boundary layer separates on the backside of the jet. Using various Eulerian and Lagrangian diagnostics, it is demonstrated that several near-wall vortical structures are formed as the wall boundary layer separates. The counter-rotating vortex pair, manifested by the presence of vortices aligned with the jet trajectory, is initiated closer to the jet exit. Moreover tornado-like wall-normal vortices originate from the separation of spanwise vorticity in the wall boundary layer at the side of the jet and from the entrainment of streamwise wall vortices in the recirculation zone on the lee side. These tornado-like vortices are absent in the case where separation is suppressed. Tornado-like vortices merge with counter-rotating vorticity originating in the jet shear layer, significantly increasing wall-normal circulation and causing deeper jet penetration into the crossflow stream.

@article{art_14,
  author     = "Schlegel, F. and Wee, D. and Marzouk, Y. M. and Ghoniem, A. F.",
  journal    = "Journal of Fluid Mechanics",
  number     = "1",
  pages      = "461-490",
  title      = "Contributions of the wall boundary layer to the formation of the counter-rotating vortex pair in transverse jets.",
  volume     = "676",
  year       = "2011",
  doi        = "10.1017/jfm.2011.59",
  abstract   = "Using high-resolution 3-D vortex simulations, this study seeks a mechanistic understanding of vorticity dynamics in transverse jets at a finite Reynolds number. A full no-slip boundary condition, rigorously formulated in terms of vorticity generation along the channel wall, captures unsteady interactions between the wall boundary layer and the jet -- in particular, the separation of the wall boundary layer and its transport into the interior. For comparison, we also implement a reduced boundary condition that suppresses the separation of the wall boundary layer away from the jet nozzle. By contrasting results obtained with these two boundary conditions, we characterize near-field vortical structures formed as the wall boundary layer separates on the backside of the jet. Using various Eulerian and Lagrangian diagnostics, it is demonstrated that several near-wall vortical structures are formed as the wall boundary layer separates. The counter-rotating vortex pair, manifested by the presence of vortices aligned with the jet trajectory, is initiated closer to the jet exit. Moreover tornado-like wall-normal vortices originate from the separation of spanwise vorticity in the wall boundary layer at the side of the jet and from the entrainment of streamwise wall vortices in the recirculation zone on the lee side. These tornado-like vortices are absent in the case where separation is suppressed. Tornado-like vortices merge with counter-rotating vorticity originating in the jet shear layer, significantly increasing wall-normal circulation and causing deeper jet penetration into the crossflow stream.",
  keywords   = "jets; vortex flows"
}

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.

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

Terrorist attacks using an aerosolized pathogen have gained credibility as a national security concern after the anthrax attacks of 2001. Inferring some important details of the attack quickly, for example, the number of people infected, the time of infection, and a representative dose received can be crucial to planning a medical response. We use a Bayesian approach, based on a short time series of diagnosed patients, to estimate a joint probability density for these parameters. We first test the formulation with idealized cases and then apply it to realistic scenarios, including the Sverdlovsk anthrax outbreak of 1979. We also use simulated outbreaks to explore the impact of model error, as when the model used for generating simulated epidemic curves does not match the model subsequently used to characterize the attack. We find that in all cases except for the smallest attacks (fewer than 100 infected people), 3–5 days of data are sufficient to characterize the outbreak to a specificity that is useful for directing an emergency response.

@article{art_12,
  author     = "Ray, J. and Marzouk, Y. M. and Najm, H. N.",
  journal    = "Statistics in Medicine",
  number     = "2",
  pages      = "101-126",
  title      = "A Bayesian approach for estimating bioterror attacks from patient data",
  volume     = "30",
  year       = "2010",
  doi        = "10.1002/sim.4090",
  abstract   = "Terrorist attacks using an aerosolized pathogen have gained credibility as a national security concern after the anthrax attacks of 2001. Inferring some important details of the attack quickly, for example, the number of people infected, the time of infection, and a representative dose received can be crucial to planning a medical response. We use a Bayesian approach, based on a short time series of diagnosed patients, to estimate a joint probability density for these parameters. We first test the formulation with idealized cases and then apply it to realistic scenarios, including the Sverdlovsk anthrax outbreak of 1979. We also use simulated outbreaks to explore the impact of model error, as when the model used for generating simulated epidemic curves does not match the model subsequently used to characterize the attack. We find that in all cases except for the smallest attacks (fewer than 100 infected people), 3--5 days of data are sufficient to characterize the outbreak to a specificity that is useful for directing an emergency response.",
  keywords   = "Bayesian inference;anthrax;Sverdlovsk outbreak;bioterrorism"
}

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.

@article{art_11,
  author     = "Wee, D. and Marzouk, Y. M. and Schlegel, F. and Ghoniem, A. F.",
  journal    = "SIAM Journal on Scientific Computing",
  number     = "4",
  pages      = "2510-2527",
  title      = "Convergence characteristics and computational cost of two algebraic kernels in vortex methods with a tree-code algorithm",
  volume     = "31",
  year       = "2009",
  doi        = "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."
}

We consider a Bayesian approach to nonlinear inverse problems in which the unknown quantity is a spatial or temporal field, endowed with a hierarchical Gaussian process prior. Computational challenges in this construction arise from the need for repeated evaluations of the forward model (e.g., in the context of Markov chain Monte Carlo) and are compounded by high dimensionality of the posterior. We address these challenges by introducing truncated Karhunen–Loève expansions, based on the prior distribution, to efficiently parameterize the unknown field and to specify a stochastic forward problem whose solution captures that of the deterministic forward model over the support of the prior. We seek a solution of this problem using Galerkin projection on a polynomial chaos basis, and use the solution to construct a reduced-dimensionality surrogate posterior density that is inexpensive to evaluate. We demonstrate the formulation on a transient diffusion equation with prescribed source terms, inferring the spatially-varying diffusivity of the medium from limited and noisy data.

@article{art_10,
  author     = "Marzouk, Y. M. and Najm, H. N.",
  journal    = "Journal of Computational Physics",
  number     = "6",
  pages      = "1862-1902",
  title      = "Dimensionality reduction and polynomial chaos acceleration of Bayesian inference in inverse problems",
  volume     = "228",
  year       = "2009",
  doi        = "10.1016/j.jcp.2008.11.024",
  abstract   = "We consider a Bayesian approach to nonlinear inverse problems in which the unknown quantity is a spatial or temporal field, endowed with a hierarchical Gaussian process prior. Computational challenges in this construction arise from the need for repeated evaluations of the forward model (e.g., in the context of Markov chain Monte Carlo) and are compounded by high dimensionality of the posterior. We address these challenges by introducing truncated Karhunen--Loève expansions, based on the prior distribution, to efficiently parameterize the unknown field and to specify a stochastic forward problem whose solution captures that of the deterministic forward model over the support of the prior. We seek a solution of this problem using Galerkin projection on a polynomial chaos basis, and use the solution to construct a reduced-dimensionality surrogate posterior density that is inexpensive to evaluate. We demonstrate the formulation on a transient diffusion equation with prescribed source terms, inferring the spatially-varying diffusivity of the medium from limited and noisy data.",
  keywords   = "Inverse problems; Bayesian inference; Dimensionality reduction; Polynomial chaos; Markov chain Monte Carlo; Galerkin projection; Gaussian processes; Karhunen--Loève expansion; RKHS"
}

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.

@article{art_9,
  author     = "Marzouk, Y. M. and Xiu, D.",
  journal    = "Communications in Computational Physics",
  number     = "1",
  pages      = "826-847",
  title      = "A stochastic collocation approach to Bayesian inference in inverse problems",
  volume     = "6",
  year       = "2009",
  doi        = "DOI: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."
}

We demonstrate the use of multiwavelet spectral polynomial chaos techniques for uncertainty quantification in non-isothermal ignition of a methane–air system. We employ Bayesian inference for identifying the probabilistic representation of the uncertain parameters and propagate this uncertainty through the ignition process. We analyze the time evolution of moments and probability density functions of the solution. We also examine the role and significance of dependence among the uncertain parameters. We finish with a discussion of the role of non-linearity and the performance of the algorithm.

@article{art_8,
  author     = "Najm, H. N. and Debusschere, B. J. and Marzouk, Y. M. and Widmer, S. and LeMaître, O.",
  journal    = "International Journal for Numerical Methods in Engineering",
  number     = "6-7",
  pages      = "789-814",
  title      = "Uncertainty quantification in chemical systems",
  volume     = "80",
  year       = "2009",
  doi        = "10.1002/nme.2551",
  abstract   = "We demonstrate the use of multiwavelet spectral polynomial chaos techniques for uncertainty quantification in non-isothermal ignition of a methane--air system. We employ Bayesian inference for identifying the probabilistic representation of the uncertain parameters and propagate this uncertainty through the ignition process. We analyze the time evolution of moments and probability density functions of the solution. We also examine the role and significance of dependence among the uncertain parameters. We finish with a discussion of the role of non-linearity and the performance of the algorithm.",
  keywords   = "uncertainty quantification; polynomial chaos; multiwavelets; chemistry; ignition"
}

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.

@article{art_7,
  author     = "Sargsyan, K. and Debusschere, B. J. and Najm, H. N. and Marzouk, Y. M.",
  journal    = "Journal of Computational and Theoretical Nanoscience",
  number     = "10",
  pages      = "2283-2297",
  title      = "Bayesian inference of spectral expansions for predictability assessment in stochastic reaction networks",
  volume     = "6",
  year       = "2009",
  doi        = "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.",
  keywords   = "uncertainty quantification; bayesian inference; polynomial chaos; stochastic reaction networks; domain decomposition; predictability"
}

We present a reformulation of the Bayesian approach to inverse problems, that seeks to accelerate Bayesian inference by using polynomial chaos (PC) expansions to represent random variables. Evaluation of integrals over the unknown parameter space is recast, more efficiently, as Monte Carlo sampling of the random variables underlying the PC expansion. We evaluate the utility of this technique on a transient diffusion problem arising in contaminant source inversion. The accuracy of posterior estimates is examined with respect to the order of the PC representation, the choice of PC basis, and the decomposition of the support of the prior. The computational cost of the new scheme shows significant gains over direct sampling.

@article{art_6,
  author     = "Marzouk, Y. M. and Najm, H. N. and Rahn, L. A.",
  journal    = "Journal of Computational Physics",
  number     = "2",
  pages      = "560-586",
  title      = "Stochastic spectral methods for efficient Bayesian solution of inverse problems",
  volume     = "224",
  year       = "2007",
  doi        = "10.1016/j.jcp.2006.10.010",
  abstract   = "We present a reformulation of the Bayesian approach to inverse problems, that seeks to accelerate Bayesian inference by using polynomial chaos (PC) expansions to represent random variables. Evaluation of integrals over the unknown parameter space is recast, more efficiently, as Monte Carlo sampling of the random variables underlying the PC expansion. We evaluate the utility of this technique on a transient diffusion problem arising in contaminant source inversion. The accuracy of posterior estimates is examined with respect to the order of the PC representation, the choice of PC basis, and the decomposition of the support of the prior. The computational cost of the new scheme shows significant gains over direct sampling.",
  keywords   = "Inverse problems; Bayesian inference; Polynomial chaos; Monte Carlo; Markov chain Monte Carlo; Spectral methods; Galerkin projection; Diffusive transport"
}

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.

@article{art_5,
  author     = "Marzouk, Y. M. and Ghoniem, A. F.",
  journal    = "Journal of Fluid Mechanics",
  number     = "1",
  pages      = "267-305",
  title      = "Vorticity structure and evolution in a transverse jet",
  volume     = "575",
  year       = "2007",
  doi        = "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."
}

A number of complex physical problems can be approached through N-body simulation, from fluid flow at high Reynolds number to gravitational astrophysics and molecular dynamics. In all these applications, direct summation is prohibitively expensive for large N and thus hierarchical methods are employed for fast summation. This work introduces new algorithms, based on k-means clustering, for partitioning parallel hierarchical N-body interactions. We demonstrate that the number of particle–cluster interactions and the order at which they are performed are directly affected by partition geometry. Weighted k-means partitions minimize the sum of clusters’ second moments and create well-localized domains, and thus reduce the computational cost of N-body approximations by enabling the use of lower-order approximations and fewer cells. We also introduce compatible techniques for dynamic load balancing, including adaptive scaling of cluster volumes and adaptive redistribution of cluster centroids. We demonstrate the performance of these algorithms by constructing a parallel treecode for vortex particle simulations, based on the serial variable-order Cartesian code developed by Lindsay and Krasny [Journal of Computational Physics 172 (2) (2001) 879–907]. The method is applied to vortex simulations of a transverse jet. Results show outstanding parallel efficiencies even at high concurrencies, with velocity evaluation errors maintained at or below their serial values; on a realistic distribution of 1.2 million vortex particles, we observe a parallel efficiency of 98% on 1024 processors. Excellent load balance is achieved even in the face of several obstacles, such as an irregular, time-evolving particle distribution containing a range of length scales and the continual introduction of new vortex particles throughout the domain. Moreover, results suggest that k-means yields a more efficient partition of the domain than a global oct-tree.

@article{art_4,
  author     = "Marzouk, Y. M. and Ghoniem, A. F.",
  journal    = "Journal of Computational Physics",
  number     = "2",
  pages      = "493-528",
  title      = "K-means clustering for optimal partitioning and dynamic load balancing of parallel hierarchical N-body simulations",
  volume     = "207",
  year       = "2005",
  doi        = "10.1016/j.jcp.2005.01.021",
  abstract   = "A number of complex physical problems can be approached through N-body simulation, from fluid flow at high Reynolds number to gravitational astrophysics and molecular dynamics. In all these applications, direct summation is prohibitively expensive for large N and thus hierarchical methods are employed for fast summation. This work introduces new algorithms, based on k-means clustering, for partitioning parallel hierarchical N-body interactions. We demonstrate that the number of particle--cluster interactions and the order at which they are performed are directly affected by partition geometry. Weighted k-means partitions minimize the sum of clusters’ second moments and create well-localized domains, and thus reduce the computational cost of N-body approximations by enabling the use of lower-order approximations and fewer cells. We also introduce compatible techniques for dynamic load balancing, including adaptive scaling of cluster volumes and adaptive redistribution of cluster centroids. We demonstrate the performance of these algorithms by constructing a parallel treecode for vortex particle simulations, based on the serial variable-order Cartesian code developed by Lindsay and Krasny [Journal of Computational Physics 172 (2) (2001) 879--907]. The method is applied to vortex simulations of a transverse jet. Results show outstanding parallel efficiencies even at high concurrencies, with velocity evaluation errors maintained at or below their serial values; on a realistic distribution of 1.2 million vortex particles, we observe a parallel efficiency of 98% on 1024 processors. Excellent load balance is achieved even in the face of several obstacles, such as an irregular, time-evolving particle distribution containing a range of length scales and the continual introduction of new vortex particles throughout the domain. Moreover, results suggest that k-means yields a more efficient partition of the domain than a global oct-tree.",
  keywords   = "k-means clustering; Treecode; N-body problems; Hierarchical methods; Parallel processing; Load balancing; Particle methods; Vortex methods; Three-dimensional flow; Transverse jet"
}

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.

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

Premixed flames encounter gradients of mixture equivalence ratio in stratified charge engines, lean premixed gas-turbine engines, and a variety of other applications. In cases for which the scales—spatial or temporal—of fuel concentration gradients in the reactants are comparable to flame scales, changes in burning rate, flammability limits, and flame structure have been observed. This paper uses an unsteady strained flame in the stagnation point configuration to examine the effect of temporal gradients on combustion in a premixed methane/air mixture. An inexact Newton backtracking method, coupled with a preconditioned Krylov subspace iterative solver, was used to improve the efficiency of the numerical solution and expand its domain of convergence in the presence of detailed chemistry. Results indicate that equivalence ratio variations with timescales lower than 10 ms have significant effects on the burning process, including reaction zone broadening, burning rate enhancement, and extension of the flammability limit toward learner mixtures. While the temperature of a flame processing a stoichiometric-to-lean equivalence ratio gradient decreased slightly within the front side of the reaction zone, radical concentrations remained elevated over the entire flame structure. These characteristics are linked to a feature reminiscent of “back-supported” flames—flames in which a stream of products resulting from burning at higher equivalence ratio is continuously supplied to lower equivalence ratio reactants. The relevant feature is the establishment of a positive temperature gradient on the products side of the flame which maintains the temperature high enough and the radical concentration sufficient to sustain combustion there. Unsteadiness in equivalence ratio produces similar gradients within the flame structure, thus compensating for the change in temperature at the leading edge of the reaction zone and accounting for an observed “flame inertia”. For sufficiently large equivalence ratio gradients, a flame starting in a stoichiometric mixture can burn through a very lean one by taking advantage of this mechanism.

@article{art_2,
  author     = "Marzouk, Y. M. and Ghoniem, A. F. and Najm, H. N.",
  journal    = "Proceedings of the Combustion Institute",
  number     = "2",
  pages      = "1859-1866",
  title      = "Dynamic response of strained premixed flames to equivalence ratio gradients",
  volume     = "28",
  year       = "2000",
  doi        = "10.1016/S0082-0784(00)80589-5",
  abstract   = "Premixed flames encounter gradients of mixture equivalence ratio in stratified charge engines, lean premixed gas-turbine engines, and a variety of other applications. In cases for which the scales—spatial or temporal—of fuel concentration gradients in the reactants are comparable to flame scales, changes in burning rate, flammability limits, and flame structure have been observed. This paper uses an unsteady strained flame in the stagnation point configuration to examine the effect of temporal gradients on combustion in a premixed methane/air mixture. An inexact Newton backtracking method, coupled with a preconditioned Krylov subspace iterative solver, was used to improve the efficiency of the numerical solution and expand its domain of convergence in the presence of detailed chemistry. Results indicate that equivalence ratio variations with timescales lower than 10 ms have significant effects on the burning process, including reaction zone broadening, burning rate enhancement, and extension of the flammability limit toward learner mixtures. While the temperature of a flame processing a stoichiometric-to-lean equivalence ratio gradient decreased slightly within the front side of the reaction zone, radical concentrations remained elevated over the entire flame structure. These characteristics are linked to a feature reminiscent of “back-supported” flames—flames in which a stream of products resulting from burning at higher equivalence ratio is continuously supplied to lower equivalence ratio reactants. The relevant feature is the establishment of a positive temperature gradient on the products side of the flame which maintains the temperature high enough and the radical concentration sufficient to sustain combustion there. Unsteadiness in equivalence ratio produces similar gradients within the flame structure, thus compensating for the change in temperature at the leading edge of the reaction zone and accounting for an observed “flame inertia”. For sufficiently large equivalence ratio gradients, a flame starting in a stoichiometric mixture can burn through a very lean one by taking advantage of this mechanism."
}

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.

@article{art_1,
  author     = "Marzouk, Y. M. and Hart, D. P.",
  journal    = "Experiments in Fluids",
  number     = "5-6",
  pages      = "401-408",
  title      = "Asymmetric autocorrelation function to resolve directional ambiguity in PIV images",
  volume     = "25",
  year       = "1998",
  doi        = "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."
}