Machine learning promises to accelerate the pace of engineering with surrogate models that mimic traditional simulation at a fraction of the cost. However, estimating the surrogate’s error remains a challenge; time-to-solution is irrelevant if the solution is too wrong. We develop two frames for understanding this challenge: (i) characterizing surrogate modeling as curve fitting and (ii) distinguishing rational and empirical error estimation approaches. The frames yield insights: accurate approximation requires sufficient resolution; rational error estimates are impractical; and empirical error estimates are uncertain. We apply these frames to three common surrogates: projection-based reduced-order models, Gaussian process emulators, and neural networks. We conclude with ideas for increasing empirical rigor in surrogate modeling.

We introduce the Lipschitz matrix: a generalization of the scalar Lipschitz constant for functions with many inputs. Among the Lipschitz matrices compatible a particular function, we choose the smallest such matrix in the Frobenius norm to encode the structure of this function. The Lipschitz matrix then provides a function-dependent metric on the input space. Altering this metric to reflect a particular function improves the performance of many tasks in computational science. Compared to the Lipschitz constant, the Lipschitz matrix reduces the worst-case cost of approximation, integration, and optimization; if the Lipschitz matrix is low-rank, this cost no longer depends on the dimension of the input, but instead on the rank of the Lipschitz matrix defeating the curse of dimensionality. Both the Lipschitz constant and matrix define uncertainty away from point queries of the function and by using the Lipschitz matrix we can reduce uncertainty. If we build a minimax space-filling design of experiments in the Lipschitz matrix metric, we can further reduce this uncertainty. When the Lipschitz matrix is approximately low-rank, we can perform parameter reduction by constructing a ridge approximation whose active subspace is the span of the dominant eigenvectors of the Lipschitz matrix. In summary, the Lipschitz matrix provides a new tool for analyzing and performing parameter reduction in complex models arising in computational science.

Parameter reduction can enable otherwise infeasible design and uncertainty studies with modern computational science models that contain several input parameters. In statistical regression, techniques for sufficient dimension reduction (SDR) use data to reduce the predictor dimension of a regression problem. A computational scientist hoping to use SDR for parameter reduction encounters a problem: a computer prediction is best represented by a deterministic function of the inputs, so data comprised of computer simulation queries fail to satisfy the SDR assumptions. To address this problem, we interpret SDR methods sliced inverse regression (SIR) and sliced average variance estimation (SAVE) as estimating the directions of a ridge function, which is a composition of a low-dimensional linear transformation with a nonlinear function. Within this interpretation, SIR and SAVE estimate matrices of integrals whose column spaces are contained in the ridge directions' span; we analyze and numerically verify convergence of these column spaces as the number of computer model queries increases. Moreover, we show example functions that are not ridge functions but whose inverse conditional moment matrices are low-rank. Consequently, the computational scientist should beware when using SIR and SAVE for parameter reduction, since SIR and SAVE may mistakenly suggest that truly important directions are unimportant.

We perform a Bayesian calibration of the freestream velocity and density starting from measurements of the pressure and heat flux at the stagnation point of a hypersonic high-enthalpy flow around a cylinder. The objective is to explore the possibility of using stagnation heat flux measurements, together with pressure measurements, to rebuild freestream conditions since such measurements are available for recent space missions but not exploited for freestream characterization. First, we formulate an algorithm of mesh adaptation, enabling accurate numerical solutions in an automatic way for a given set of inputs. Secondly, active subspaces are used to find a low-dimensional dependence structures in the input-to-output map of the forward numerical solver. Then, surrogate models on the active variables are used to accelerate the forward uncertainty propagation by Monte Carlo sampling and the Markov Chain Monte Carlo sampling of the posterior distribution for Bayesian inversion. A preliminary sensitivity analysis with sparse Polynomial Dimensional Decomposition is performed on the chemical model of the air mixture, to determine the most influential uncertain chemical parameters in the forward problem. Then, the forward and backward methodologies are applied to the simulation of a hypersonic flow around a cylinder, in conditions for which experimental data are available, revealing new insights towards the potential exploitation of heat flux data for freestream rebuilding.

Multivariate functions encountered in high-dimensional uncertainty quantification problems often vary most strongly 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 the case 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.

Assuming a multivariate function is Lipschitz continuous is one way to arrive at the curse of dimensionality---the situation where the cost of tasks such as approximation, integration, and optimization grow exponentially in the input dimension of the function. Here we introduce the Lipschitz matrix, a generalization of the scalar Lipschitz constant that alters the norm in which the distance between points is measured. When this Lipschitz matrix is low-rank, we are able to avoid this curse by noting that the quantity of interest is a ridge function with a low-dimensional subspace (the range of the Lipschitz matrix) onto which we can restrict the input without loss. Even if the Lipschitz matrix is not low rank, it can still reduce the constant associated with these tasks compared the Lipschitz constant when the singular values of the Lipschitz matrix decay. In addition, the Lipschitz matrix has other uses: it can identify an active subspace for parameter space dimension reduction; it can define uncertainty and provide informative bounds in high-dimensional spaces; and it motivates a space filling design of experiments to minimize uncertainty. As there are a limited number of situations where the Lipschitz matrix can be identified analytically, we show that the Lipschitz matrix can be estimated by solving a semidefinite program using a finite number of function samples or gradients.

Consider an open set $\mathbb{D}\subseteq\mathbb{R}^n$, equipped with a probability measure $\mu$. An important characteristic of a smooth function $f:\mathbb{D}\rightarrow\mathbb{R}$ is its \emph{second-moment matrix} $\Sigma_{\mu}:=\int \nabla f(x) \nabla f(x)^* \mu(dx) \in\mathbb{R}^{n\times n}$, where $\nabla f(x)\in\mathbb{R}^n$ is the gradient of $f(\cdot)$ at $x\in\mathbb{D}$ and $*$ stands for transpose. For instance, the span of the leading $r$ eigenvectors of $\Sigma_{\mu}$ forms an \emph{active subspace} of $f(\cdot)$, which contains the directions along which $f(\cdot)$ changes the most and is of particular interest in \emph{ridge approximation}. In this work, we propose a simple algorithm for estimating $\Sigma_{\mu}$ from random point evaluations of $f(\cdot)$ \emph{without} imposing any structural assumptions on $\Sigma_{\mu}$. Theoretical guarantees for this algorithm are established with the aid of the same technical tools that have proved valuable in the context of covariance matrix estimation from partial measurements.

Sufficient dimension reduction (SDR) provides a framework for reducing the predictor space dimension in statistical regression problems. We consider SDR in the context of dimension reduction for deterministic functions of several variables such as those arising in computer experiments. In this context, SDR can reveal low-dimensional ridge structure in functions. Two algorithms for SDR—sliced inverse regression (SIR) and sliced average variance estimation (SAVE)—approximate matrices of integrals using a sliced mapping of the response. We interpret this sliced approach as a Riemann sum approximation of the particular integrals arising in each algorithm. We employ the well-known tools from numerical analysis—namely, multivariate numerical integration and orthogonal polynomials—to produce new algorithms that improve upon the Riemann sum-based numerical integration in SIR and SAVE. We call the new algorithms Lanczos–Stieltjes inverse regression (LSIR) and Lanczos–Stieltjes average variance estimation (LSAVE) due to their connection with Stieltjes’ method—and Lanczos’ related discretization—for generating a sequence of polynomials that are orthogonal with respect to a given measure. We show that this approach approximates the desired integrals, and we study the behavior of LSIR and LSAVE with two numerical examples. The quadrature-based LSIR and LSAVE eliminate the first-order algebraic convergence rate bottleneck resulting from the Riemann sum approximation, thus enabling high-order numerical approximations of the integrals when appropriate. Moreover, LSIR and LSAVE perform as well as the best-case SIR and SAVE implementations (e.g., adaptive partitioning of the response space) when low-order numerical integration methods (e.g., simple Monte Carlo) are used.

Many of the input-parameter-to-output-quantity-of-interest maps that arise in computational science admit a surprising low-dimensional structure, where the outputs vary primarily along a handful of directions in the high-dimensional input space. This type of structure is well modeled by a ridge function, which is a composition of a low-dimensional linear transformation with a nonlinear function. If the goal is to compute statistics of the output-e.g., as in uncertainty quantification or robust design-then one should exploit this low-dimensional structure, when present, to accelerate computations. We develop Gaussian quadrature and the associated polynomial approximation for one-dimensional ridge functions. The key elements of our method are (i) approximating the univariate density of the given linear combination of inputs by repeated convolutions and (ii) a Lanczos Stieltjes method for constructing orthogonal polynomials and Gaussian quadrature.

A modified, deterministic SEIR model is developed for the 2014 Ebola epidemic occurring in the West African nations of Guinea, Liberia, and Sierra Leone. The model describes the dynamical interaction of susceptible and infected populations, while accounting for the effects of hospitalization and the spread of disease through interactions with deceased, but infectious, individuals. Using data from the World Health Organization (WHO), parameters within the model are fit to recent estimates of infected and deceased cases from each nation. The model is then analyzed using these parameter values. Finally, several metrics are proposed to determine which of these nations is in greatest need of additional resources to combat the spread of infection. These include local and global sensitivity metrics of both the infected population and the basic reproduction number with respect to rates of hospitalization and proper burial.

Knowledge of accurate values of elastic modulus of (Al 1-x Sc x )N is required for design of piezoelectric resonators and related devices.Thin films of (Al 1-x Sc x )N across the entire composition space are deposited and characterized.Accuracy of modulus measurements is improved and quantified by removing the influence of substrate effects and by direct comparison of experimental results with density functional theory calculations.The 5%-30% Sc compositional range is of particular interest for piezoelectric applications and is covered at higher compositional resolution here than in previous work.The reduced elastic modulus is found to decrease by as much as 40% with increasing Sc concentration in the wurtzite phase according to both experimental and computational techniques, whereas Sc-rich rocksalt-structured films exhibit little variation in modulus with composition.

Inexpensive surrogates are useful for reducing the cost of science and engineering studies involving large-scale, complex computational models with many input parameters. A ridge approximation is one class of surrogate that models a quantity of interest as a nonlinear function of a few linear combinations of the input parameters. When used in parameter studies (e.g., optimization or uncertainty quantification), ridge approximations allow the low-dimensional structure to be exploited, reducing the effective dimension. We introduce a new, fast algorithm for constructing a ridge approximation where the nonlinear function is a polynomial. This polynomial ridge approximation is chosen to minimize least squares mismatch between the surrogate and the quantity of interest on a given set of inputs. Naively, this would require optimizing both the polynomial coefficients and the linear combination of weights, the latter of which define a low-dimensional subspace of the input space. However, given a fixed subspace the optimal polynomial can be found by solving a linear leastsquares problem. Hence using variable projection the polynomial can be implicitly defined, leaving an optimization problem over the subspace alone. Here we develop an algorithm that finds this polynomial ridge approximation by minimizing over the Grassmann manifold of low-dimensional subspaces using a Gauss{Newton method. Our Gauss{Newton method has superior theoretical guarantees and faster convergence on our numerical examples than the alternating approach for polynomial ridge approximation earlier proposed by Constantine, Eftekhari, Hokanson, and Ward [Comput. Methods Appl. Mech. Engrg., 326 (2017), pp. 402{421] that alternates between (i) optimizing the polynomial coefficients given the subspace and (ii) optimizing the subspace given the coefficients.

Many of the input-parameter-to-output-quantity-of-interest maps that arise in computational science admit a surprising low-dimensional structure, where the outputs vary primarily along a handful of directions in the high-dimensional input space. This type of structure is well modeled by a ridge function, which is a composition of a low-dimensional linear transformation with a nonlinear function. If the goal is to compute statistics of the output (e.g., as in uncertainty quantification or robust design) then one should exploit this low-dimensional structure, when present, to accelerate computations. We develop Gaussian quadrature and the associated polynomial approximation for one-dimensional ridge functions. The key elements of our method are (i) approximating the univariate density of the given linear combination of inputs by repeated convolutions and (ii) a Lanczos-Stieltjes method for constructing orthogonal polynomials and Gaussian quadrature.

Problems in uncertainty quantification (UQ) suffer from the curse of dimensionality. One approach to address this issue is to identify and exploit low-dimensional structure in the underlying model. In this paper, we show that low-dimensional active subspaces are not present in a particular model arising from acoustic scattering. This suggests that UQ techniques based on active subspaces are not appropriate for this problem.

Abstract. Active subspaces can effectively reduce the dimension of high-dimensional parameter studies enabling otherwise infeasible experiments with expensive simulations. The key components of active subspace methods are the eigenvectors of a symmetric, positive semidefinite matrix whose elements are the average products of partial derivatives of the simulation’s input/output map. We study a Monte Carlo method for approximating the eigenpairs of this matrix. We offer both theoretical results based on recent non-asymptotic random matrix theory and a practical approach based on the bootstrap. We extend the analysis to the case when the gradients are approximated, for example, with finite differences. Our goal is to provide guidance for two questions that arise in active subspaces: (i) How many gradient samples does one need to accurately approximate the eigenvalues and subspaces? (ii) What can be said about the accuracy of the estimated subspace, both theoretically and practically? We test the approach on both simple quadratic functions where the active subspace is known and a parameterized PDE with 100 variables characterizing the coefficients of the differential operator.

Design and optimization benefit from understanding the dependence of a quantity of interest (e.g., a design objective or constraint function) on the design variables. A low-dimensional active subspace, when present, identifies important directions in the space of design variables; perturbing a design along the active subspace associated with a particular quantity of interest changes that quantity more, on average, than perturbing the design orthogonally to the active subspace. This low-dimensional structure provides insights that characterize the dependence of quantities of interest on design variables. Airfoil design in a transonic flow field with a parameterized geometry is a popular test problem for design methodologies. We examine two particular airfoil shape parameterizations, PARSEC and CST, and study the active subspaces present in two common design quantities of interest, transonic lift and drag coefficients, under each shape parameterization. We mathematically relate the two parameterizations with a common polynomial series. The active subspaces enable low-dimensional approximations of lift and drag that relate to physical airfoil properties. In particular, we obtain and interpret a two-dimensional approximation of both transonic lift and drag, and we show how these approximation inform a multi-objective design problem.

We propose a gradient-based method for detecting and exploiting low-dimensional input parameter dependence of multivariate functions. The methodology consists in minimizing an upper bound, obtained by Poincar\'e-type inequalities, on the approximation error. The resulting method can be used to approximate vector-valued functions (e.g., functions taking values in $\mathbb{R}^n$ or functions taking values in function spaces) and generalizes the notion of active subspaces associated with scalar-valued functions. A comparison with the truncated Karhunen-Lo\`eve decomposition shows that using gradients of the function can yield more effective dimension reduction. Numerical examples reveal that the choice of norm on the codomain of the function can have a significant impact on the function's low-dimensional approximation.

The hydrology of high-elevation, mountainous regions is poorly represented in Earth Systems Models (ESMs), yet these ecosystems play an important role in the storage and land-atmosphere exchange of water. As much of the western United States' water comes from water stored in the snowpack (snow water equivalent, SWE), model representation of these regions is important. This study assesses how uncertainty in both model parameters and forcing affect simulated snow processes through sensitivity analysis (active subspaces) on model inputs (meteorological forcing and model input parameters) for a widely used snow model. Observations from an AmeriFlux tower at the Niwot Ridge research site are used to force an integrated, single-column hydrologic model, ParFlow-CLM. This study finds that trees can mute the effects of snow albedo causing the evergreen needleleaf scenarios to be sensitive primarily to hydrologic forcing while bare ground simulations are more sensitive to the snow parameters. The bare ground scenarios are most sensitive overall. Both forcing and model input parameters are important for obtaining accurate hydrologic model results.