Paul G. Constantine

Paul G. Constantine, Michael S. Eldred, Eric T. Phipps

Multivariate global polynomial approximations - such as polynomial chaos or stochastic collocation methods - are now in widespread use for sensitivity analysis and uncertainty quantification. The pseudospectral variety of these methods uses a numerical integration rule to approximate the Fourier-type coefficients of a truncated expansion in orthogonal polynomials. For problems in more than two or three dimensions, a sparse grid numerical integration rule offers accuracy with a smaller node set compared to tensor product approximation. However, when using a sparse rule to approximately integrate these coefficients, one often finds unacceptable errors in the coefficients associated with higher degree polynomials. By reexamining Smolyak's algorithm and exploiting the connections between interpolation and projection in tensor product spaces, we construct a sparse pseudospectral approximation method that accurately reproduces the coefficients of basis functions that naturally correspond to the sparse grid integration rule. The compelling numerical results show that this is the proper way to use sparse grid integration rules for pseudospectral approximation.

Andrew T. Glaws, Paul G. Constantine, R. Dennis Cook

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.

Paul G. Constantine, Carson Kent, Tan Bui-Thanh

The Markov chain Monte Carlo (MCMC) method is the computational workhorse for Bayesian inverse problems. However, MCMC struggles in high-dimensional parameter spaces, since its iterates must sequentially explore the high-dimensional space. This struggle is compounded in physical applications when the nonlinear forward model is computationally expensive. One approach to accelerate MCMC is to reduce the dimension of the state space. Active subspaces are part of an emerging set of tools for subspace-based dimension reduction. An active subspace in a given inverse problem indicates a separation between a low-dimensional subspace that is informed by the data and its orthogonal complement that is constrained by the prior. With this information, one can run the sequential MCMC on the active variables while sampling independently according to the prior on the inactive variables. However, this approach to increase efficiency may introduce bias. We provide a bound on the Hellinger distance between the true posterior and its active subspace- exploiting approximation. And we demonstrate the active subspace-accelerated MCMC on two computational examples: (i) a two-dimensional parameter space with a quadratic forward model and one-dimensional active subspace and (ii) a 100-dimensional parameter space with a PDE-based forward model and a two-dimensional active subspace.

Paul G. Constantine, Alireza Doostan

Renewable energy researchers use computer simulation to aid the design of lithium ion storage devices. The underlying models contain several physical input parameters that affect model predictions. Effective design and analysis must understand the sensitivity of model predictions to changes in model parameters, but global sensitivity analyses become increasingly challenging as the number of input parameters increases. Active subspaces are part of an emerging set of tools for discovering and exploiting low-dimensional structures in the map from high-dimensional inputs to model outputs. We extend linear and quadratic model-based heuristic for active sub- space discovery to time-dependent processes and apply the resulting technique to a lithium ion battery model. The results reveal low-dimensional structure and sensitivity metrics that a designer may exploit to study the relationship between parameters and predictions.

Paul G. Constantine, Zachary del Rosario, Gianluca Iaccarino

A ridge function is a function of several variables that is constant along certain directions in its domain. Using classical dimensional analysis, we show that many physical laws are ridge functions; this fact yields insight into the structure of physical laws and motivates further study into ridge functions and their properties. We also connect dimensional analysis to modern subspace-based techniques for dimension reduction, including active subspaces in deterministic approximation and sufficient dimension reduction in statistical regression.

Zachary J. Grey, Paul G. Constantine

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.

Jeffrey M. Hokanson, Paul G. Constantine

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 squared 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 least-squares problem, and hence by using variable projection the polynomial can be implicitly found leaving an optimization problem over the subspace alone. We provide an algorithm that finds this polynomial ridge approximation by minimizing over the Grassmann manifold of low-dimensional subspaces using a Gauss-Newton method. We provide details of this optimization algorithm and demonstrate its performance on several numerical examples. Our Gauss-Newton method has superior theoretical guarantees and faster convergence than the alternating approach for polynomial ridge approximation earlier proposed by Constantine, Eftekhari, Hokanson, and Ward [https://doi.org/10.1016/j.cma.2017.07.038] that alternates between (i) optimizing the polynomial coefficients given the subspace and (ii) optimizing the subspace given the coefficients.

Paul G. Constantine, Brian Zaharatos, Mark Campanelli

Predictions from science and engineering models depend on the values of the model's input parameters. As the number of parameters increases, algorithmic parameter studies like optimization or uncertainty quantification require many more model evaluations. One way to combat this curse of dimensionality is to seek an alternative parameterization with fewer variables that produces comparable predictions. The active subspace is a low-dimensional linear subspace defined by important directions in the model's input space; input perturbations along these directions change the model's prediction more, on average, than perturbations orthogonal to the important directions. We describe a method for checking if a model admits an exploitable active subspace, and we apply this method to a single-diode solar cell model with five input parameters. We find that the maximum power of the solar cell has a dominant one-dimensional active subspace, which enables us to perform thorough parameter studies in one dimension instead of five.

Paul G. Constantine, Qiqi Wang

We present a method for approximating the solution of a parameterized, nonlinear dynamical system using an affine combination of solutions computed at other points in the input parameter space. The coefficients of the affine combination are computed with a nonlinear least squares procedure that minimizes the residual of the governing equations. The approximation properties of this residual minimizing scheme are comparable to existing reduced basis and POD-Galerkin model reduction methods, but its implementation requires only independent evaluations of the nonlinear forcing function. It is particularly appropriate when one wishes to approximate the states at a few points in time without time marching from the initial conditions. We prove some interesting characteristics of the scheme including an interpolatory property, and we present heuristics for mitigating the effects of the ill-conditioning and reducing the overall cost of the method. We apply the method to representative numerical examples from kinetics - a three state system with one parameter controlling the stiffness - and conductive heat transfer - a nonlinear parabolic PDE with a random field model for the thermal conductivity.

Paul G. Constantine, Eric T. Phipps

We seek to approximate a composite function h(x) = g(f(x)) with a global polynomial. The standard approach chooses points x in the domain of f and computes h(x) at each point, which requires an evaluation of f and an evaluation of g. We present a Lanczos-based procedure that implicitly approximates g with a polynomial of f. By constructing a quadrature rule for the density function of f, we can approximate h(x) using many fewer evaluations of g. The savings is particularly dramatic when g is much more expensive than f or the dimension of x is large. We demonstrate this procedure with two numerical examples: (i) an exponential function composed with a rational function and (ii) a Navier-Stokes model of fluid flow with a scalar input parameter that depends on multiple physical quantities.

Andrew Glaws, Paul G. Constantine, R. Dennis Cook

We investigate the application of sufficient dimension reduction (SDR) to a noiseless data set derived from a deterministic function of several variables. In this context, SDR provides a framework for ridge recovery. In this second part, we explore the numerical subtleties associated with using two inverse regression methods---sliced inverse regression (SIR) and sliced average variance estimation (SAVE)---for ridge recovery. This includes a detailed numerical analysis of the eigenvalues of the resulting matrices and the subspaces spanned by their columns. After this analysis, we demonstrate the methods on several numerical test problems.

Andrew Glaws, Paul G. Constantine

Sufficient dimension reduction (SDR) provides a framework for reducing the predictor space dimension in regression problems. We consider SDR in the context of deterministic functions of several variables such as those arising in computer experiments. In this context, SDR serves as a methodology for uncovering ridge structure in functions, and two primary 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 well-known tools from numerical analysis---namely, multivariate tensor product Gauss-Christoffel quadrature 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 to a given measure. We show that the quadrature-based approach approximates the desired integrals, and we study the behavior of LSIR and LSAVE with three numerical examples. As expected in high order numerical integration, the quadrature-based LSIR and LSAVE exhibit exponential convergence in the integral approximations compared to the first order convergence of the classical SIR and SAVE. The disadvantage of LSIR and LSAVE is that the underlying tensor product quadrature suffers from the curse of dimensionality---that is, the number of quadrature nodes grows exponentially with the input space dimension. Therefore, the proposed approach is most appropriate for deterministic functions with fewer than ten independent inputs.

Paul G. Constantine, Qiqi Wang

This manuscript is superseded by Constantine, Dow, and Wang's "Active Subspaces in Theory and Practice: Applications to Kriging Surfaces" [SIAM J. of Sci. Comput., 36 (2014), pp. A1500-A1524]. Many multivariate functions encountered in practice vary primarily along a few directions in the space of input parameters. When these directions correspond with coordinate directions, one may apply global sensitivity measures to determine the parameters with the greatest contribution to the function's variability. However, these methods perform poorly when the directions of variability are not aligned with the natural coordinates of the input space. We present a method for detecting the directions of variability of a function using evaluations of its derivative with respect to the input parameters. We demonstrate how to exploit these directions to construct a surrogate function that depends on fewer variables than the original function, thus reducing the dimension of the original problem. We apply this procedure to an exercise in uncertainty quantification using an elliptic PDE with a model for the coefficients that depends on 250 independent parameters. The dimension reduction procedure identifies a 5-dimensional subspace suitable for constructing surrogates.

Paul G. Constantine, David F. Gleich, Gianluca Iaccarino

Recent work has explored solver strategies for the linear system of equations arising from a spectral Galerkin approximation of the solution of PDEs with parameterized (or stochastic) inputs. We consider the related problem of a matrix equation whose matrix and right hand side depend on a set of parameters (e.g. a PDE with stochastic inputs semidiscretized in space) and examine the linear system arising from a similar Galerkin approximation of the solution. We derive a useful factorization of this system of equations, which yields bounds on the eigenvalues, clues to preconditioning, and a flexible implementation method for a wide array of problems. We complement this analysis with (i) a numerical study of preconditioners on a standard elliptic PDE test problem and (ii) a fluids application using existing CFD codes; the MATLAB codes used in the numerical studies are available online.

Andrew Glaws, Paul G. Constantine

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.

Paul G. Constantine, Zachary del Rosario, Gianluca Iaccarino

Classical dimensional analysis has two limitations: (i) the computed dimensionless groups are not unique, and (ii) the analysis does not measure relative importance of the dimensionless groups. We propose two algorithms for estimating unique and relevant dimensionless groups assuming the experimenter can control the system's independent variables and evaluate the corresponding dependent variable; e.g., computer experiments provide such a setting. The first algorithm is based on a response surface constructed from a set of experiments. The second algorithm uses many experiments to estimate finite differences over a range of the independent variables. Both algorithms are semi-empirical because they use experimental data to complement the dimensional analysis. We derive the algorithms by combining classical semi-empirical modeling with active subspaces, which---given a probability density on the independent variables---yield unique and relevant dimensionless groups. The connection between active subspaces and dimensional analysis also reveals that all empirical models are ridge functions, which are functions that are constant along low-dimensional subspaces in its domain. We demonstrate the proposed algorithms on the well-studied example of viscous pipe flow---both turbulent and laminar cases. The results include a new set of two dimensionless groups for turbulent pipe flow that are ordered by relevance to the system; the precise notion of relevance is closely tied to the derivative based global sensitivity metric from Sobol' and Kucherenko.

Paul G. Constantine, Armin Eftekhari, Jeffrey Hokanson et al.

Response surfaces are common surrogates for expensive computer simulations in engineering analysis. However, the cost of fitting an accurate response surface increases exponentially as the number of model inputs increases, which leaves response surface construction intractable for high-dimensional, nonlinear models. We describe ridge approximation for fitting response surfaces in several variables. A ridge function is constant along several directions in its domain, so fitting occurs on the coordinates of a low-dimensional subspace of the input space. We review essential theory for ridge approximation---e.g., the best mean-squared approximation and an optimal low-dimensional subspace---and we prove that the gradient-based active subspace is near-stationary for the least-squares problem that defines an optimal subspace. Motivated by the theory, we propose a computational heuristic that uses an estimated active subspace as an initial guess for a ridge approximation fitting problem. We show a simple example where the heuristic fails, which reveals a type of function for which the proposed approach is inappropriate. We then propose a simple alternating heuristic for fitting a ridge function, and we demonstrate the effectiveness of the active subspace initial guess applied to an airfoil model of drag as a function of its 18 shape parameters.

Andrew Glaws, Paul G. Constantine, John Shadid et al.

Magnetohydrodynamics (MHD)---the study of electrically conducting fluids---can be harnessed to produce efficient, low-emissions power generation. Today, computational modeling assists engineers in studying candidate designs for such generators. However, these models are computationally expensive, so studying the effects of the model's many input parameters on output predictions is typically infeasible. We study two approaches for reducing the input dimension of the models: (i) classical dimensional analysis based on the inputs' units and (ii) active subspaces, which reveal low-dimensional subspaces in the space of inputs that affect the outputs the most. We also review the mathematical connection between the two approaches that leads to consistent application. The dimension reduction yields insights into the driving factors in the MHD power generation models. We study both the simplified Hartmann problem, which admits closed form expressions for the quantities of interest, and a large-scale computational model with adjoint capabilities that enable the derivative computations needed to estimate the active subspaces.

Paul G. Constantine, Paul Diaz

Predictions from science and engineering models depend on several input parameters. Global sensitivity analysis quantifies the importance of each input parameter, which can lead to insight into the model and reduced computational cost; commonly used sensitivity metrics include Sobol' total sensitivity indices and derivative-based global sensitivity measures. Active subspaces are an emerging set of tools for identifying important directions in a model's input parameter space; these directions can be exploited to reduce the model's dimension enabling otherwise infeasible parameter studies. In this paper, we develop global sensitivity metrics called activity scores from the active subspace, which yield insight into the important model parameters. We mathematically relate the activity scores to established sensitivity metrics, and we discuss computational methods to estimate the activity scores. We show two numerical examples with algebraic functions taken from simplified engineering models. For each model, we analyze the active subspace and discuss how to exploit the low-dimensional structure. We then show that input rankings produced by the activity scores are consistent with rankings produced by the standard metrics.

Paul G. Constantine, Armin Eftekhari, Michael B. Wakin

Active subspaces are an emerging set of tools for identifying and exploiting the most important directions in the space of a computer simulation's input parameters; these directions depend on the simulation's quantity of interest, which we treat as a function from inputs to outputs. To identify a function's active subspace, one must compute the eigenpairs of a matrix derived from the function's gradient, which presents challenges when the gradient is not available as a subroutine. We numerically study two methods for estimating the necessary eigenpairs using only linear measurements of the function's gradient. In practice, these measurements can be estimated by finite differences using only two function evaluations, regardless of the dimension of the function's input space.

Paul G. Constantine, Eric Dow, Qiqi Wang

Many multivariate functions in engineering models vary primarily along a few directions in the space of input parameters. When these directions correspond to coordinate directions, one may apply global sensitivity measures to determine the most influential parameters. However, these methods perform poorly when the directions of variability are not aligned with the natural coordinates of the input space. We present a method to first detect the directions of the strongest variability using evaluations of the gradient and subsequently exploit these directions to construct a response surface on a low-dimensional subspace---i.e., the active subspace---of the inputs. We develop a theoretical framework with error bounds, and we link the theoretical quantities to the parameters of a kriging response surface on the active subspace. We apply the method to an elliptic PDE model with coefficients parameterized by 100 Gaussian random variables and compare it with a local sensitivity analysis method for dimension reduction.

Paul G. Constantine, David F. Gleich, Gianluca Iaccarino

We apply polynomial approximation methods -- known in the numerical PDEs context as spectral methods -- to approximate the vector-valued function that satisfies a linear system of equations where the matrix and the right hand side depend on a parameter. We derive both an interpolatory pseudospectral method and a residual-minimizing Galerkin method, and we show how each can be interpreted as solving a truncated infinite system of equations; the difference between the two methods lies in where the truncation occurs. Using classical theory, we derive asymptotic error estimates related to the region of analyticity of the solution, and we present a practical residual error estimate. We verify the results with two numerical examples.