Search

This work develops an efficient real-time inverse formulation for inferring the aerodynamic surface pressures on a hypersonic vehicle from sparse measurements of the structural strain. The approach aims to provide real-time estimates of the aerodynamic loads acting on the vehicle for ground and flight testing, as well as guidance, navigation, and control applications. Specifically, the approach targets hypersonic flight conditions where direct measurement of the surface pressures is challenging due to the harsh aerothermal environment. For problems employing a linear elastic structural model, we show that the inference problem can be posed as a least-squares problem with a linear constraint arising from a finite element discretization of the governing elasticity partial differential equation. Due to the linearity of the problem, an explicit solution is given by the normal equations. Pre-computation of the resulting inverse map enables rapid evaluation of the surface pressure and corresponding integrated quantities, such as the force and moment coefficients. The inverse approach additionally allows for uncertainty quantification, providing insights for theoretical recoverability and robustness to sensor noise. Numerical studies demonstrate the estimator performance for reconstructing the surface pressure field, as well as the force and moment coefficients, for the Initial Concept 3.X (IC3X) conceptual hypersonic vehicle., Comment: 22 pages, 15 figures

This work considers the computation of risk measures for quantities of interest governed by PDEs with Gaussian random field parameters using Taylor approximations. While efficient, Taylor approximations are local to the point of expansion, and hence may degrade in accuracy when the variances of the input parameters are large. To address this challenge, we approximate the underlying Gaussian measure by a mixture of Gaussians with reduced variance in a dominant direction of parameter space. Taylor approximations are constructed at the means of each Gaussian mixture component, which are then combined to approximate the risk measures. The formulation is presented in the setting of infinite-dimensional Gaussian random parameters for risk measures including the mean, variance, and conditional value-at-risk. We also provide detailed analysis of the approximations errors arising from two sources: the Gaussian mixture approximation and the Taylor approximations. Numerical experiments are conducted for a semilinear advection-diffusion-reaction equation with a random diffusion coefficient field and for the Helmholtz equation with a random wave speed field. For these examples, the proposed approximation strategy can achieve less than $1\%$ relative error in estimating CVaR with only $\mathcal{O}(10)$ state PDE solves, which is comparable to a standard Monte Carlo estimate with $\mathcal{O}(10^4)$ samples, thus achieving significant reduction in computational cost. The proposed method can therefore serve as a way to rapidly and accurately estimate risk measures under limited computational budgets., Comment: 34 Pages, 13 Figures, 1 Table

This paper establishes optimal convergence rates for estimation of structured covariance operators of Gaussian processes. We study banded operators with kernels that decay rapidly off-the-diagonal and $L^q$-sparse operators with an unordered sparsity pattern. For these classes of operators, we find the minimax optimal rate of estimation in operator norm, identifying the fundamental dimension-free quantities that determine the sample complexity. In addition, we prove that tapering and thresholding estimators attain the optimal rate. The proof of the upper bound for tapering estimators requires novel techniques to circumvent the issue that discretization of a banded operator does not result, in general, in a banded covariance matrix. To derive lower bounds for banded and $L^q$-sparse classes, we introduce a general framework to lift theory from high-dimensional matrix estimation to the operator setting. Our work contributes to the growing literature on operator estimation and learning, building on ideas from high-dimensional statistics while also addressing new challenges that emerge in infinite dimension., Comment: 38 pages, 3 figures

We present a scalable and efficient framework for the inference of spatially-varying parameters of continuum materials from image observations of their deformations. Our goal is the nondestructive identification of arbitrary damage, defects, anomalies and inclusions without knowledge of their morphology or strength. Since these effects cannot be directly observed, we pose their identification as an inverse problem. Our approach builds on integrated digital image correlation (IDIC, Besnard Hild, Roux, 2006), which poses the image registration and material inference as a monolithic inverse problem, thereby enforcing physical consistency of the image registration using the governing PDE. Existing work on IDIC has focused on low-dimensional parameterizations of materials. In order to accommodate the inference of heterogeneous material propertes that are formally infinite dimensional, we present $\infty$-IDIC, a general formulation of the PDE-constrained coupled image registration and inversion posed directly in the function space setting. This leads to several mathematical and algorithmic challenges arising from the ill-posedness and high dimensionality of the inverse problem. To address ill-posedness, we consider various regularization schemes, namely $H^1$ and total variation for the inference of smooth and sharp features, respectively. To address the computational costs associated with the discretized problem, we use an efficient inexact-Newton CG framework for solving the regularized inverse problem. In numerical experiments, we demonstrate the ability of $\infty$-IDIC to characterize complex, spatially varying Lam\'e parameter fields of linear elastic and hyperelastic materials. Our method exhibits (i) the ability to recover fine-scale and sharp material features, (ii) mesh-independent convergence performance and hyperparameter selection, (iii) robustness to observational noise.

We present an efficient and scalable algorithm for performing matrix-vector multiplications ("matvecs") for block Toeplitz matrices. Such matrices, which are shift-invariant with respect to their blocks, arise in the context of solving inverse problems governed by autonomous systems, and time-invariant systems in particular. In this article, we consider inverse problems that are solved for inferring unknown parameters from observational data of a linear time-invariant dynamical system given in the form of partial differential equations (PDEs). Matrix-free Newton-conjugate-gradient methods are often the gold standard for solving these inverse problems, but they require numerous actions of the Hessian on a vector. Matrix-free adjoint-based Hessian matvecs require solution of a pair of linearized forward/adjoint PDE solves per Hessian action, which may be prohibitive for large-scale inverse problems, especially when efficient low-rank approximations of the Hessian are not readily available, such as for hyperbolic PDE operators. Time invariance of the forward PDE problem leads to a block Toeplitz structure of the discretized parameter-to-observable (p2o) map defining the mapping from inputs (parameters) to outputs (observables) of the PDEs. This block Toeplitz structure enables us to exploit two key properties: (1) compact storage of the p2o map and its adjoint; and (2) efficient fast Fourier transform (FFT)-based Hessian matvecs. The proposed algorithm is mapped onto large multi-GPU clusters and achieves more than 80 percent of peak bandwidth on an NVIDIA A100 GPU. Excellent weak scaling is shown for up to 48 A100 GPUs. For the targeted problems, the implementation executes Hessian matvecs within fractions of a second, orders of magnitude faster than can be achieved by the conventional matrix-free Hessian matvecs via forward/adjoint PDE solves.

Estimating the parameters of compact binaries which coalesce and produce gravitational waves is a challenging Bayesian inverse problem. Gravitational-wave parameter estimation lies within the class of multifidelity problems, where a variety of models with differing assumptions, levels of fidelity, and computational cost are available for use in inference. In an effort to accelerate the solution of a Bayesian inverse problem, cheaper surrogates for the best models may be used to reduce the cost of likelihood evaluations when sampling the posterior. Importance sampling can then be used to reweight these samples to represent the true target posterior, incurring a reduction in the effective sample size. In cases when the problem is high dimensional, or when the surrogate model produces a poor approximation of the true posterior, this reduction in effective samples can be dramatic and render multifidelity importance sampling ineffective. We propose a novel method of tempered multifidelity importance sampling in order to remedy this issue. With this method the biasing distribution produced by the low-fidelity model is tempered, allowing for potentially better overlap with the target distribution. There is an optimal temperature which maximizes the efficiency in this setting, and we propose a low-cost strategy for approximating this optimal temperature using samples from the untempered distribution. In this paper, we motivate this method by applying it to Gaussian target and biasing distributions. Finally, we apply it to a series of problems in gravitational wave parameter estimation and demonstrate improved efficiencies when applying the method to real gravitational wave detections., Comment: 19 pages, 10 figures, 1 table

This paper studies sparse covariance operator estimation for nonstationary Gaussian processes with sharply varying marginal variance and small correlation lengthscale. We introduce a covariance operator estimator that adaptively thresholds the sample covariance function using an estimate of the variance components. Building on recent results from empirical process theory, we derive an operator norm bound on the estimation error in terms of the sparsity level of the covariance and the expected supremum of the normalized process. Our theory and numerical simulations demonstrate the advantage of adaptive threshold estimators over universal threshold and sample covariance estimators in nonstationary settings., Comment: 38 pages, 5 figures

This is a concise mathematical introduction to Monte Carlo methods, a rich family of algorithms with far-reaching applications in science and engineering. Monte Carlo methods are an exciting subject for mathematical statisticians and computational and applied mathematicians: the design and analysis of modern algorithms are rooted in a broad mathematical toolbox that includes ergodic theory of Markov chains, Hamiltonian dynamical systems, transport maps, stochastic differential equations, information theory, optimization, Riemannian geometry, and gradient flows, among many others. These lecture notes celebrate the breadth of mathematical ideas that have led to tangible advancements in Monte Carlo methods and their applications. To accommodate a diverse audience, the level of mathematical rigor varies from chapter to chapter, giving only an intuitive treatment to the most technically demanding subjects. The aim is not to be comprehensive or encyclopedic, but rather to illustrate some key principles in the design and analysis of Monte Carlo methods through a carefully-crafted choice of topics that emphasizes timeless over timely ideas. Algorithms are presented in a way that is conducive to conceptual understanding and mathematical analysis -- clarity and intuition are favored over state-of-the-art implementations that are harder to comprehend or rely on ad-hoc heuristics. To help readers navigate the expansive landscape of Monte Carlo methods, each algorithm is accompanied by a summary of its pros and cons, and by a discussion of the type of problems for which they are most useful. The presentation is self-contained, and therefore adequate for self-guided learning or as a teaching resource. Each chapter contains a section with bibliographic remarks that will be useful for those interested in conducting research on Monte Carlo methods and their applications., Comment: 150 pages, 21 figures

We propose an operator learning approach to accelerate geometric Markov chain Monte Carlo (MCMC) for solving infinite-dimensional Bayesian inverse problems (BIPs). While geometric MCMC employs high-quality proposals that adapt to posterior local geometry, it requires repeated computations of gradients and Hessians of the log-likelihood, which becomes prohibitive when the parameter-to-observable (PtO) map is defined through expensive-to-solve parametric partial differential equations (PDEs). We consider a delayed-acceptance geometric MCMC method driven by a neural operator surrogate of the PtO map, where the proposal exploits fast surrogate predictions of the log-likelihood and, simultaneously, its gradient and Hessian. To achieve a substantial speedup, the surrogate must accurately approximate the PtO map and its Jacobian, which often demands a prohibitively large number of PtO map samples via conventional operator learning methods. In this work, we present an extension of derivative-informed operator learning [O'Leary-Roseberry et al., J. Comput. Phys., 496 (2024)] that uses joint samples of the PtO map and its Jacobian. This leads to derivative-informed neural operator (DINO) surrogates that accurately predict the observables and posterior local geometry at a significantly lower training cost than conventional methods. Cost and error analysis for reduced basis DINO surrogates are provided. Numerical studies demonstrate that DINO-driven MCMC generates effective posterior samples 3--9 times faster than geometric MCMC and 60--97 times faster than prior geometry-based MCMC. Furthermore, the training cost of DINO surrogates breaks even compared to geometric MCMC after just 10--25 effective posterior samples., Comment: Updated manuscript: changed title, changed format, typo correction, and minor terminology changes

This paper investigates covariance operator estimation via thresholding. For Gaussian random fields with approximately sparse covariance operators, we establish non-asymptotic bounds on the estimation error in terms of the sparsity level of the covariance and the expected supremum of the field. We prove that thresholded estimators enjoy an exponential improvement in sample complexity compared with the standard sample covariance estimator if the field has a small correlation lengthscale. As an application of the theory, we study thresholded estimation of covariance operators within ensemble Kalman filters., Comment: 29 pages, 2 figures

Filtering is concerned with online estimation of the state of a dynamical system from partial and noisy observations. In applications where the state of the system is high dimensional, ensemble Kalman filters are often the method of choice. These algorithms rely on an ensemble of interacting particles to sequentially estimate the state as new observations become available. Despite the practical success of ensemble Kalman filters, theoretical understanding is hindered by the intricate dependence structure of the interacting particles. This paper investigates ensemble Kalman filters that incorporate an additional resampling step to break the dependency between particles. The new algorithm is amenable to a theoretical analysis that extends and improves upon those available for filters without resampling, while also performing well in numerical examples., Comment: 32 pages, 5 figures

We present an efficient matrix-free point spread function (PSF) method for approximating operators that have locally supported non-negative integral kernels. The method computes impulse responses at scattered points, and interpolates these impulse responses to approximate integral kernel entries. Impulse responses are computed by applying the operator to Dirac comb batches of point sources, which are chosen via an ellipsoid packing procedure. Evaluation of kernel entries allows us to construct a hierarchical matrix approximation of the operator, which is used for further matrix computations. We illustrate the end-to-end method on a blur problem, then use the method to build preconditioners for the Hessian in two inverse problems governed by partial differential equations (PDEs): inversion for the basal friction coefficient in an ice sheet flow problem and for the initial condition in an advective-diffusive transport problem. While for many ill-posed inverse problems the Hessian of the data misfit term exhibits a low rank structure, and hence a low rank approximation is suitable, for many problems of practical interest the numerical rank of the Hessian is still large. But Hessian impulse responses typically become more local as the numerical rank increases, which benefits the PSF method. Numerical results reveal that the PSF preconditioner clusters the spectrum of the preconditioned Hessian near one, yielding roughly 5x-10x reductions in the required number of PDE solves, as compared to regularization preconditioning and no preconditioning. We also present a numerical study for the influence of various parameters (that control the shape of the impulse responses) on the effectiveness of the advection-diffusion Hessian approximation. The results show that the PSF-based preconditioners are able to form good approximations of high-rank Hessians using a small number of operator applications.

Identifying parameters of computational models from experimental data, or model calibration, is fundamental for assessing and improving the predictability and reliability of computer simulations. In this work, we propose a method for Bayesian calibration of models that predict morphological patterns of diblock copolymer (Di-BCP) thin film self-assembly while accounting for various sources of uncertainties in pattern formation and data acquisition. This method extracts the azimuthally-averaged power spectrum (AAPS) of the top-down microscopy characterization of Di-BCP thin film patterns as summary statistics for Bayesian inference of model parameters via the pseudo-marginal method. We derive the analytical and approximate form of a conditional likelihood for the AAPS of image data. We demonstrate that AAPS-based image data reduction retains the mutual information, particularly on important length scales, between image data and model parameters while being relatively agnostic to the aleatoric uncertainties associated with the random long-range disorder of Di-BCP patterns. Additionally, we propose a phase-informed prior distribution for Bayesian model calibration. Furthermore, reducing image data to AAPS enables us to efficiently build surrogate models to accelerate the proposed Bayesian model calibration procedure. We present the formulation and training of two multi-layer perceptrons for approximating the parameter-to-spectrum map, which enables fast integrated likelihood evaluations. We validate the proposed Bayesian model calibration method through numerical examples, for which the neural network surrogate delivers a fivefold reduction of the number of model simulations performed for a single calibration task., Comment: Minor changes from the original submission, including a change in the title

We propose a novel machine learning framework for solving optimization problems governed by large-scale partial differential equations (PDEs) with high-dimensional random parameters. Such optimization under uncertainty (OUU) problems may be computational prohibitive using classical methods, particularly when a large number of samples is needed to evaluate risk measures at every iteration of an optimization algorithm, where each sample requires the solution of an expensive-to-solve PDE. To address this challenge, we propose a new neural operator approximation of the PDE solution operator that has the combined merits of (1) accurate approximation of not only the map from the joint inputs of random parameters and optimization variables to the PDE state, but also its derivative with respect to the optimization variables, (2) efficient construction of the neural network using reduced basis architectures that are scalable to high-dimensional OUU problems, and (3) requiring only a limited number of training data to achieve high accuracy for both the PDE solution and the OUU solution. We refer to such neural operators as multi-input reduced basis derivative informed neural operators (MR-DINOs). We demonstrate the accuracy and efficiency our approach through several numerical experiments, i.e. the risk-averse control of a semilinear elliptic PDE and the steady state Navier--Stokes equations in two and three spatial dimensions, each involving random field inputs. Across the examples, MR-DINOs offer $10^{3}$--$10^{7} \times$ reductions in execution time, and are able to produce OUU solutions of comparable accuracies to those from standard PDE based solutions while being over $10 \times$ more cost-efficient after factoring in the cost of construction.

We explore using neural operators, or neural network representations of nonlinear maps between function spaces, to accelerate infinite-dimensional Bayesian inverse problems (BIPs) with models governed by nonlinear parametric partial differential equations (PDEs). Neural operators have gained significant attention in recent years for their ability to approximate the parameter-to-solution maps defined by PDEs using as training data solutions of PDEs at a limited number of parameter samples. The computational cost of BIPs can be drastically reduced if the large number of PDE solves required for posterior characterization are replaced with evaluations of trained neural operators. However, reducing error in the resulting BIP solutions via reducing the approximation error of the neural operators in training can be challenging and unreliable. We provide an a priori error bound result that implies certain BIPs can be ill-conditioned to the approximation error of neural operators, thus leading to inaccessible accuracy requirements in training. To reliably deploy neural operators in BIPs, we consider a strategy for enhancing the performance of neural operators, which is to correct the prediction of a trained neural operator by solving a linear variational problem based on the PDE residual. We show that a trained neural operator with error correction can achieve a quadratic reduction of its approximation error, all while retaining substantial computational speedups of posterior sampling when models are governed by highly nonlinear PDEs. The strategy is applied to two numerical examples of BIPs based on a nonlinear reaction--diffusion problem and deformation of hyperelastic materials. We demonstrate that posterior representations of the two BIPs produced using trained neural operators are greatly and consistently enhanced by error correction.

In this paper we propose a variant of enriched Galerkin methods for second order elliptic equations with over-penalization of interior jump terms. The bilinear form with interior over-penalization gives a non-standard norm which is different from the discrete energy norm in the classical discontinuous Galerkin methods. Nonetheless we prove that optimal a priori error estimates with the standard discrete energy norm can be obtained by combining a priori and a posteriori error analysis techniques. We also show that the interior over-penalization is advantageous for constructing preconditioners robust to mesh refinement by analyzing spectral equivalence of bilinear forms. Numerical results are included to illustrate the convergence and preconditioning results., Comment: preconditioning for anisotropic cases is improved

Many modern algorithms for inverse problems and data assimilation rely on ensemble Kalman updates to blend prior predictions with observed data. Ensemble Kalman methods often perform well with a small ensemble size, which is essential in applications where generating each particle is costly. This paper develops a non-asymptotic analysis of ensemble Kalman updates that rigorously explains why a small ensemble size suffices if the prior covariance has moderate effective dimension due to fast spectrum decay or approximate sparsity. We present our theory in a unified framework, comparing several implementations of ensemble Kalman updates that use perturbed observations, square root filtering, and localization. As part of our analysis, we develop new dimension-free covariance estimation bounds for approximately sparse matrices that may be of independent interest.

Directed self-assembly (DSA) of block-copolymers (BCPs) is one of the most promising developments in the cost-effective production of nanoscale devices. The process makes use of the natural tendency for BCP mixtures to form nanoscale structures upon phase separation. The phase separation can be directed through the use of chemically patterned substrates to promote the formation of morphologies that are essential to the production of semiconductor devices. Moreover, the design of substrate pattern can formulated as an optimization problem for which we seek optimal substrate designs that effectively produce given target morphologies. In this paper, we adopt a phase field model given by a nonlocal Cahn--Hilliard partial differential equation (PDE) based on the minimization of the Ohta--Kawasaki free energy, and present an efficient PDE-constrained optimization framework for the optimal design problem. The design variables are the locations of circular- or strip-shaped guiding posts that are used to model the substrate chemical pattern. To solve the ensuing optimization problem, we propose a variant of an inexact Newton conjugate gradient algorithm tailored to this problem. We demonstrate the effectiveness of our computational strategy on numerical examples that span a range of target morphologies. Owing to our second-order optimizer and fast state solver, the numerical results demonstrate five orders of magnitude reduction in computational cost over previous work. The efficiency of our framework and the fast convergence of our optimization algorithm enable us to rapidly solve the optimal design problem in not only two, but also three spatial dimensions., Comment: 35 Pages, 17 Figures

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

We propose derivative-informed neural operators (DINOs), a general family of neural networks to approximate operators as infinite-dimensional mappings from input function spaces to output function spaces or quantities of interest. After discretizations both inputs and outputs are high-dimensional. We aim to approximate not only the operators with improved accuracy but also their derivatives (Jacobians) with respect to the input function-valued parameter to empower derivative-based algorithms in many applications, e.g., Bayesian inverse problems, optimization under parameter uncertainty, and optimal experimental design. The major difficulties include the computational cost of generating derivative training data and the high dimensionality of the problem leading to large training cost. To address these challenges, we exploit the intrinsic low-dimensionality of the derivatives and develop algorithms for compressing derivative information and efficiently imposing it in neural operator training yielding derivative-informed neural operators. We demonstrate that these advances can significantly reduce the costs of both data generation and training for large classes of problems (e.g., nonlinear steady state parametric PDE maps), making the costs marginal or comparable to the costs without using derivatives, and in particular independent of the discretization dimension of the input and output functions. Moreover, we show that the proposed DINO achieves significantly higher accuracy than neural operators trained without derivative information, for both function approximation and derivative approximation (e.g., Gauss-Newton Hessian), especially when the training data are limited.