## COMPUTING ACTIVE SUBSPACES WITH MONTE CARLO

PAUL G. CONSTANTINE ∗ AND DAVID F. GLEICH †

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.

Key words. active subspaces, dimension reduction

## AMS subject classifications.

1. Introduction. Engineering models typically contain several input parameters that must be specified to produce a set of model outputs that contains one or more quantities of interest. The engineer's goal is to characterize the behavior of the quantities of interest as functions of the model's inputs. However, parameter studiessuch as optimization and uncertainty quantification-are challenging when the number of inputs is large and the model involves an expensive computer simulation. In such cases, the engineer may analyze the output's sensitivity with respect to inputs to identify a subset of inputs whose variation changes the outputs the most [26]. In the best case, she can then limit parameter studies to key parameters and thus reduce the dimension of the parameter study. This approach is appropriate when varying important parameters changes the outputs much more than varying the unimportant parameters. However, a model's output may depend on all the parameters through certain linear combinations, which generalizes seeking key parameters to seeking key directions in the parameter space. The active subspace identifies important directions in the model's input space with respect to a particular quantity of interest; perturbing the inputs along these important directions changes the quantity of interest more, on average, than perturbing the inputs in orthogonal directions [7]. For parameter studies whose work depends exponentially on the number of parameters-e.g., integration or response surface construction-the active subspace-enabled dimension reduction can permit otherwise infeasible studies.

The active subspace is defined by a set of eigenvectors corresponding to large eigenvalues of the average outer product of the gradient with itself. These eigenpairs are properties of the map between model inputs and outputs, like Fourier coefficients or the Lipschitz constant. To determine if a function admits a low-dimensional active subspace-and thus reduce the dimension of the parameter studies-one must estimate these eigenpairs. This estimation is problematic because the elements of the matrix defining the eigenpairs are themselves high-dimensional integrals. Most de-

∗ Department of Applied Mathematics and Statistics, Colorado School of Mines, Golden, CO 80401 ( paul.constantine@mines.edu ).

† Department of Computer Science, Purdue University, West Lafayette, IN ( dgleich@purdue.edu )

terministic numerical integration rules are impractical beyond a handful of variables, especially if the integrand is costly to evaluate. We therefore focus on a Monte Carlo approach to approximate the eigenpairs, where we take advantage of recent theoretical results that bound the number of samples needed to approximate the spectrum of sums of random matrices. Monte Carlo is attractive because it makes few restrictions on the function defining the quantity of interest. Under additional assumptions, one may be able to outperform Monte Carlo with specialized integration rules for integrands that depend on many variables, e.g., with sparse grids [5] or quasi-Monte Carlo [6].

In what follows, we analyze a Monte Carlo method for estimating the eigenpairs that uses independent samples of the function's gradient. After formally defining the active subspace in Section 2, we employ results from Tropp [30] and Gittens and Tropp [14] to bound the probability that the estimated eigenvalues deviate from the true eigenvalues, which yields lower bounds on the number of samples needed for accurate estimation. We extend these results to the case where samples are approximate gradients (e.g., finite difference approximations). In Section 4 we discuss a practical bootstrap approach to study the variability in the estimated eigenvalues, and we demonstrate these procedures numerically in Section 5.

Notation. We use bold lower case letters to denote vectors and bold upper case letters to denote matrices. Finite sample estimates are denoted with hats, e.g., ˆ C ≈ C . The functional λ k ( ) · denotes the k th eigenvalue of its argument, ordered from algebraically largest to smallest; all matrices are symmetric, so the ordering is meaningful. A λ on its own is an eigenvalue. Norms of vectors and matrices are 2-norms; the matrix 2-norm is the operator-induced norm given by the largest singular value. The partial ordering operator /precedesequal is defined as follows: A /precedesequal B means that B -A is positive semidefinite.

2. Active subspaces. We represent the map from simulation inputs to the scalar-valued quantity of interest by a function f : X → R , where X ⊆ R m , with m &gt; 1, represents the set of simulation inputs, which we assume is centered at the origin and scaled so that each component of x ∈ X has the same range. Let R m be equipped with a weight function ρ : R m → R + that is bounded, strictly positive on the domain X , and zero outside of X . We also assume that ρ is both separable and normalized to integrate to 1. In the context of uncertainty quantification, this weight function represents a given probability density function on the inputs; examples in this context include Gaussian, uniform, or data-conditioned Bayesian posterior density functions. We assume f is differentiable and absolutely continuous, and we denote the gradient ∇ x f ( x ) = [ ∂f/∂x , . . . , ∂f/∂x 1 m T ] oriented as a column vector.

We are interested in the following matrix, denoted C and defined as

$$C \, = \, \int ( \nabla _ { x } f ) ( \n$$

Samarov studied this matrix as one of several average derivative functionals in the context of regression, where f is the regression function [27]. The matrix C is symmetric and positive semi-definite, so it has a real eigenvalue decomposition

$$\begin{array} {$$

where W is the orthogonal matrix of eigenvectors. Partition the eigenpairs,

$$W = \left [ W _ { 1 } \ \ W _ { 2 } \right ], \$$

where W 1 contains the first n &lt; m eigenvectors, and Λ 1 contains the n largest eigenvalues. The eigenvectors define new coordinates

$$y = W _ { 1 } ^ { T } x, \quad z = W _ { 2 }$$

We call the column space of W 1 the active subspace and the corresponding y the active variables . Similarly, W 2 defines the inactive subspace with corresponding inactive variables z . The following two lemmas justify this characterization; these are proved in [7].

Lemma 2.1. The mean-squared directional derivative of f with respect to the eigenvector w i is equal to the corresponding eigenvalue,

$$\int \left ( ( \nabla _ { \mathbf x } f ) ^ { T$$

Lemma 2.2. The mean-squared gradients of f with respect to the coordinates y and z satisfy

$$\dots \dots \dots _ { \lambda \beta } \nu \\ \int ( \nabla _ { y } f ) ^ { T } ( \nabla _ { y } f ) \, \rho \, d x = & \, \text{trace} \left ( W _ { 1 } ^ { T } C W _ { 1 } \right ) = \lambda _ { 1 } + \dots + \lambda _ { n }, \\ \int ( \nabla _ { z } f ) ^ { T } ( \nabla _ { z } f ) \, \rho \, d x = & \, \text{trace} \left ( W _ { 2 } ^ { T } C W _ { 2 } \right ) = \lambda _ { n + 1 } + \dots + \lambda _ { m }. \\ \text{The eigenvalues $\Lambda$ and eigenvectors $W$ are properties of $f$. If the $m-n$ trailing} \\ \dots \dots \dots \lambda \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \end{cases}$$

The eigenvalues Λ and eigenvectors W are properties of f . If the m -n trailing eigenvalues Λ 2 are exactly zero, then f is constant along the directions corresponding to W 2 . If Λ 2 is not exactly zero but significantly smaller than Λ 1 , then f changes less, on average, in response to small changes in z than small changes in y . If f admits such a property, we would like to discover and exploit it in parameter studies by focusing on the variables y . In other words, we can reduce the dimension of the parameter studies from m to n &lt; m .

Two special cases illustrate the active subspace. The first class of functions are index models that have the form f ( x ) = h ( A T x ), where A ∈ R m k × and h : R k → R . In this case, C has rank at most k , and the active subspace is a subspace of the range of A . If k = 1, then the one-dimesional active subspace can be discovered with a single evaluation of ∇ x f at any x ∈ X such that the derivative h ′ ( A T x ) is not zero. The second special case is a function of the form f ( x ) = h ( x T H x ) / 2, where h : R → R , and H is a symmetric m × m matrix. In this case

$$C \, = \, H \left ( \int ( h ^ { \prime } ) ^ { 2 } \, \text{xx} ^ { T } \, \rho \, d \text{x} \right ) H ^ { T }, \, \text{quad \quad } \, ( 2. 7 ) \\ _ { T \, \text{w} \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \dots \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \dots \, \dots \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \dots \, \delta \, \dots \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \dots \, \dots \, \dots \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \dots \, \delta \, \delta \, \dots \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \dots \, \delta \, \dots \, \dots \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \dots \, \dots \, \delta \, \dots \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \dots \, \dots \, \dots \, \dots \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \dots \, \delta \, \delta \, \delta \, \dots \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \dots \, \delta \, \delta \, \dots \, \dots \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \dots \, \delta \, \dots \, \dots \, \dots \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \dots \, \dots \, \delta \, \delta \, \dots \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \dots \, \dots \, \delta \, \dots \, \dots \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \dots \, \dots \, \dots \, \delta \, \dots \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dd \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \d$$

where h ′ = h ′ ( x T H x ) is the derivative of h . This implies that the null space of C is the null space of H provided that h ′ is non-degenerate. We study the example where h t ( ) = t in Section 5.

If we can estimate Λ and W from (2.2), then we can approximate f ( x ) with a model of the form

$$f ( \mathbf x ) \, \approx \, g ( W _ { 1 } ^ { T } \mathbf x ), \quad \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \$$

where g : R n → R is an appropriately constructed map. In [7], we derive error bounds for this approximation with a particular choice of g . We extend those error bounds to the case when W 1 is estimated with some error. The main goal of this paper is to study the error in W 1 when C is estimated with Monte Carlo.

- 2.1. Related literature. The idea of studying the eigenpairs of the average outer product of the gradients arose in statistics as average derivative functionals [27, 18] for exploring structure in regression functions. In contrast to our f ( x ), the regression function is generally unknown; to estimate the gradients of the unknown regression function, one can first fit a kernel-based model to a set of predictor/response pairs and then compute gradients from the approximation [31, 12]. In our case, the function is a map between the inputs and outputs of an engineering simulation; there is no random noise as in the regression problem. The set up in Russi's Ph.D. thesis [25] is closer to ours. He applies the methods to physical simulations of chemical kinetics; this work is where we encountered the term active subspace . Recent work in approximation theory by Fornasier, Schnass, and Vybiral [11] attempts to discover the parameters of the active subspace solely through queries of the function; guarantees on reconstruction follow from compressed sensing results under the assumption that f is an index model.

If the matrix C were given as an input matrix, or if we could easily compute matrix-vector products with C , then we could employ recent procedures for randomized low-rank approximation to estimate the desired eigenpairs [17, 13]-assuming C is well approximated by a low-rank matrix, which is often the case in practice. Unfortunately, we do not have easy access to the elements of C ; estimating its eigenpairs requires estimating its elements. There may be fruitful relationships with low-rank approximation of quasimatrices and cmatrices [29] that are worth exploring.

- 3. Computing active subspaces. If drawing independent samples from the density ρ is cheap and simple, then a straightforward and easy-to-implement random sampling method to approximate the eigenvalues Λ and eigenvectors W proceeds as follows.
- 1. Draw N samples x j independently from the measure ρ .
- 2. For each x j , compute ∇ x f j = ∇ x f ( x j ).
- 3. Approximate

$$C \, & \approx \, \hat { C } \, = \, \frac { 1 } { N } \sum _ { j = 1 } ^ { N } ( \nabla _ { \mathbf x } f _ { j } ) ( \nabla _ { \mathbf x } f _ { j } ) ^ { T }. & & ( 3. 1 ) \\. \quad. \quad. \quad. \quad. \quad \hat { \wedge } \quad \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wangle } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat{ \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { \wedge } \hat { u }$$

- 4. Compute the eigendecomposition ˆ = C ˆ ˆ ˆ W W Λ T .

The last step is equivalent to computing the full SVD of the matrix

$$\hat { B } \, = \, \frac { 1 } { \sqrt { N } } \left [ \nabla _ { x } f _ { 1 } \ \cdots \ \nabla _ { x } f _ { N } \right ] \, = \, \hat { W } \hat { \Sigma } \hat { V } ^ { T }, \quad & ( 3. 2 ) \\ \text{andard manipulations show that the } \hat { \Lambda } = \hat { \Sigma } \hat { \Sigma } ^ { T }. \text{ and the left singular vectors}$$

For many simulations, the number m of input parameters is small enough (e.g., tens to thousands) that computing the full eigendecomposition (3.1) or singular value decomposition (3.2) is negligible compared to the cost of computing the gradient N times; we consider this to be our case of interest. We are therefore concerned with understanding the number of gradient samples needed so that the estimates Λ and ˆ ˆ W are close to the true Λ and W .

where standard manipulations show that the Λ = Σ ˆ ˆ ˆ Σ T , and the left singular vectors are the desired eigenvectors. The SVD perspective was developed by Russi [25] as the method to discover the active subspace. This SVD perspective also calls to mind randomized methods for subsampling the columns of ˆ , where B N /greatermuch m [13]. If it were possible to evaluate the importance of a column of ˆ B without explicitly computing ∇ x f ( x ), then such ideas might prove useful.

We apply recent work by Tropp [30] and Gittens and Tropp [14] on the spectrum of sums of random matrices to answer these questions. We were motivated to use these tools by Section 7 in Gittens and Tropp [14], which studies the spectrum of a finite sample estimate of a covariance matrix for a Gaussian random vector. In the present case, the gradient vector ∇ x f ( x ) is a deterministic function of x . However, if x j are drawn independently at random according to the density ρ , then we can interpret ∇ x f ( x j ) as a random draw from an unknown density. This is a standard interpretation of Monte Carlo techniques for integration [24]. In principle, our analysis approach could apply to model reduction of high-dimensional systems that use Grammian matrices [1].

Theorem 3.1. Assume that ‖∇ x f ‖ ≤ L for all x ∈ X . Then for 0 &lt; ε ≤ 1 ,

$$\mathbb { P } \left \{ \hat { \lambda } _ { k } \geq ( 1 + \varepsilon ) \lambda _ { k } \right \} \, \leq \, ( m - k + 1 ) \exp \left ( \frac { - N \lambda _ { k } \varepsilon ^ { 2 } } { 4 L ^ { 2 } } \right ), \quad \ \ ( 3. 3 ) \$$

and

$$\mathbb { P } \left \{ \hat { \lambda } _ { k } \leq ( 1 - \varepsilon ) \lambda _ { k } \right \} \, \leq \, k \exp \left ( \frac { - N \lambda _ { k } ^ { 2 } \varepsilon ^ { 2 } } { 4 \lambda _ { 1 } L ^ { 2 } } \right ).$$

The key to establishing Theorem 3.1 is a matrix Bernstein inequality from Theorem 5.3 in Gittens and Tropp [14]. When we apply this concentration result, we set

$$X _ { j } = \nabla _ { \mathbf x } f _ { j } \nabla _ { \mathbf x } f _ { j } ^ { T }.$$

Thus, each X j is an independent random sample of a matrix from the same distribution. Under this notion of randomness,

$$\mathbb { E } \left [ X _ { j } \right ] = \int \nabla _ { x } f _ { j } \nabla _ { x } f _ { j } ^ { T } \, \rho \, d x = \int \nabla _ { x } f \nabla _ { x } f ^ { T } \, \rho \, d x = C. \quad \quad ( 3. 6 ) \\. \quad \quad \quad$$

For completeness, we restate Theorem 5.3 from [14].

Theorem 3.2 ( Eigenvalue Bernstein Inequality for Subexponential Matrices, Theorem 5.3 [14] ). Consider a finite sequence { X j } of independent, random, selfadjoint matrices with dimension n , all of which satisfy the subexponential moment growth condition

$$\mathbb { E } \left [ X _ { j } ^ { m } \right ] \preceq \frac { m! } { 2 } B ^ { m - 2 } \Sigma _ { j } ^ { 2 } \quad \text{ for } m = 2, 3, 4, \dots \\ \colon a \text{ positive constant and } \Sigma _ { i } ^ { 2 } \text{ are positive-semidefinite matrices.}$$

where B is a positive constant and Σ 2 j are positive-semidefinite matrices. Given an integer k ≤ n , set

$$\mu _ { k } = \lambda _ { k } \left ( \sum _ { j } \mathbb { E } \left [ X _ { j } \right ] \right ). \\.$$

Choose V + as an orthogonal matrix of size n × n -k +1 that satisfies

$$\mu _ { k } = \lambda _ { \max } \left ( \sum _ { j } V _ { + } ^ { T } ( \mathbb { E } X _ { j } ) V _ { + } \right ),$$

and define

$$\sigma _ { k } ^ { 2 } = \lambda _ { \max } \left ( \sum _ { j } V _ { + } ^ { T } \Sigma _ { j } ^ { 2 } V _ { + } \right ).$$

Then, for any τ ≥ 0 ,

$$\Phi \left \{ \lambda _ { k } \left ( \sum _ { j } X _ { j } \right ) \geq \mu _ { k } + \tau \right \} \leq \begin{cases} ( n - k + 1 ) \exp ( - \tau ^ { 2 } / ( 4 \sigma _ { k } ^ { 2 } ) ), & \tau \leq \sigma _ { k } ^ { 2 } / B, \\ ( n - k + 1 ) \exp ( - \tau / ( 4 B ) ), & \tau \geq \sigma _ { k } ^ { 2 } / B. \end{cases} \\ \text{Proof. (Theorem 3.1.)} \, \text{We begin with the upper estimate (3.3). First note that}$$

Proof . (Theorem 3.1.) We begin with the upper estimate (3.3). First note that

$$\mathbb { P } \left \{ \lambda _ { k } ( \hat { C } ) \geq \lambda _ { k } ( C ) + t \right \} \, = \, \mathbb { P } \left \{ \lambda _ { k } \left ( \sum _ { j = 1 } ^ { N } \nabla _ { \mathbf x } f _ { j } \nabla _ { \mathbf x } f _ { j } ^ { T } \right ) \geq N \lambda _ { k } + N t \right \}. \quad ( 3. 7 ) \\ \text{In this form we can apply Theorem 3.2.} \, \text{We check that the bound on the gradient's}$$

  In this form we can apply Theorem 3.2. We check that the bound on the gradient's norm implies that the matrix ∇ x f ∇ x f T satisfies the subexponential growth condition:

Next we set

$$\text{orm implies that the matrix } V _ { x } J V _ { x } f ^ { * } \text{ sustines the subsequonential growth condition} \text{:} \\ \int \left ( \nabla _ { x } f \nabla _ { x } f ^ { T } \right ) ^ { p } \rho \, d x & = \int ( \nabla _ { x } f ^ { T } \nabla _ { x } f ) ^ { p - 1 } \, \nabla _ { x } f \nabla _ { x } f ^ { T } \rho \, d x \\ & \preceq \left ( L ^ { 2 } \right ) ^ { p - 1 } \int \nabla _ { x } f \nabla _ { x } f ^ { T } \rho \, d x \\ & \preceq \frac { p! } { 2 } \, \left ( L ^ { 2 } \right ) ^ { p - 2 } \, ( L ^ { 2 } C ). \\ \text{text we set} & \quad \cdot \quad \cdot$$

$$\mu _ { k } \, = \, \lambda _ { k } \left ( \sum _ { j = 1 } ^ { N } \int \nabla _ { x } f _ { j } \nabla _ { x } f _ { j } ^ { T } \, \rho \, d x \right ) \, = \, N \lambda _ { k }, \quad \, ( 3. 9 ) \\ \text{e simplified using the identically distributed samples of } x _ { j }. \text{ Choose } W _ { + } =$$

where we simplified using the identically distributed samples of x j . Choose W + = W (: , k : m ) to be the last m -k +1 eigenvectors of C , and note that

$$\lambda _ { \max } \left ( \sum _ { j = 1 } ^ { N } W _ { + } ^ { T } \left ( \int \nabla _ { x } f \nabla _ { x } f ^ { T } \rho \, d x \right ) W _ { + } \right ) \, = \, N \lambda _ { \max } ( W _ { + } ^ { T } C W _ { + } ) \, = \, N \lambda _ { k } \, = \, \mu _ { k }, \\ \text{as required by Twoorem 3.2. Define}$$

as required by Theorem 3.2. Define

$$\sigma _ { k } ^ { 2 } \, = \, \lambda _ { \max } \left ( \sum _ { j = 1 } ^ { N } W _ { + } ^ { T } ( L ^ { 2 } C ) W _ { + } \right ) \, = \, N L ^ { 2 } \, \lambda _ { \max } \left ( W _ { + } ^ { T } C W _ { + } \right ) \, = \, N L ^ { 2 } \lambda _ { k }. \quad ( 3. 1 1 ) \\ \text{With these quantities, Theorem 3.2 states}$$

With these quantities, Theorem 3.2 states

$$\mathbb { P } \left \{ \lambda _ { k } \left ( \sum _ { j = 1 } ^ { N } \nabla _ { \mathbf x } f _ { j } \nabla _ { \mathbf x } f _ { j } ^ { T } \right ) \geq N \lambda _ { k } + N t \right \} \, \leq \, ( m - k + 1 ) \exp \left ( \frac { - ( N t ) ^ { 2 } } { 4 \sigma _ { k } ^ { 2 } } \right ) \quad ( 3. 1 2 ) \\ \text{when } N t \leq \sigma _ { k } ^ { 2 } / L ^ { 2 }. \text{ Applying this theorem with } t = \varepsilon \lambda _ { k }, \, \varepsilon \leq 1, \text{ and the computed}$$

  when Nt ≤ σ /L 2 k 2 . Applying this theorem with t = ελ k , ε ≤ 1, and the computed σ 2 k = NL λ 2 k yields the upper estimate (3.3).

For the lower estimate,

$$& \text{For the lower estimate,} \\ & \mathbb { P } \left \{ \lambda _ { k } ( \tilde { C } ) \leq \lambda _ { k } ( C ) - t \right \} \\ & = \mathbb { P } \left \{ - \lambda _ { k } ( \tilde { C } ) \geq - \lambda _ { k } ( C ) + t \right \} \\ & = \mathbb { P } \left \{ - \lambda _ { k } \left ( \sum _ { j = 1 } ^ { N } \nabla _ { x } f _ { j } \nabla _ { x } f _ { j } ^ { T } \right ) \geq - N \lambda _ { k } ( C ) + N t \right \} \\ & = \mathbb { P } \left \{ \lambda _ { m - k + 1 } \left ( \sum _ { j = 1 } ^ { N } \left ( - \nabla _ { x } f _ { j } \nabla _ { x } f _ { j } ^ { T } \right ) \right ) \geq N \lambda _ { m - k + 1 } ( - C ) + N t \right \} \\ & = \mathbb { P } \left \{ \lambda _ { k ^ { \prime } } \left ( \sum _ { j = 1 } ^ { N } \left ( - \nabla _ { x } f _ { j } \nabla _ { x } f _ { j } ^ { T } \right ) \right ) \geq N \lambda _ { k ^ { \prime } } ( - C ) + N t \right \}, \\ \text{for} \, k ^ { \prime } = m - k + 1. \, \text{We can now apply Theorem 3.2 again.  The subexponential growth} \\ \text{condition is satisfied since} \\ & \int \left ( - \nabla _ { \cdot \cdot \cdot } f \nabla _ { \cdot \cdot } f ^ { T } \right \} _ { \ a d x } \, \rightarrow \, \int \left ( \nabla _ { \cdot \cdot } f \nabla _ { \cdot \cdot } f ^ { T } \right \} _ { \ a d x } \, \rightarrow \, \underline { p }! \, \left ( T ^ { 2 } \right \} ^ { p - 2 } \left ( T ^ { 2 } C \right \rangle \, \right ) \, \left ( 3 \, 1 4 \right \rangle$$

for k ′ = m -k +1. We can now apply Theorem 3.2 again. The subexponential growth condition is satisfied since

Set

$$\int \left ( - \nabla _ { x } f \nabla _ { x } f ^ { T } \right ) ^ { p } \rho \, d x \, \preceq \, \int \left ( \nabla _ { x } f \nabla _ { x } f ^ { T } \right ) ^ { p } \rho \, d x \, \preceq \, \frac { p! } { 2 } \, \left ( L ^ { 2 } \right ) ^ { p - 2 } \left ( L ^ { 2 } C \right ). \quad \left ( 3. 1 4 \right ) \\ \text{Set}$$

$$\mu _ { k ^ { \prime } } \, = \, \lambda _ { k ^ { \prime } } \left ( \sum _ { j = 1 } ^ { N } \int ( - \nabla _ { x } f _ { j } \nabla _ { x } f _ { j } ^ { T } ) \, \rho \, d x \right ) \, = \, N \lambda _ { k ^ { \prime } } ( - C ). \quad \ \ ( 3. 1 5 ) \\ \intertext { t } W _ { + } = W ( \colon, 1 \, \colon k ) \, \text{ to be the first $k$ eigenvectors of $C$, and note that}$$

Set W + = W (: , 1 : k ) to be the first k eigenvectors of C , and note that

$$\text{vew } \text{$\mathcal{r} + - \mathcal{r} \left \langle, 1 \cdot \mathcal{w} \right \rangle \, \text{or use any negative} \\ \lambda _ { \max } \left ( \sum _ { j = 1 } ^ { N } W _ { + } ^ { T } \left ( \int \left ( - \nabla _ { x } f _ { j } \nabla _ { x } f _ { j } ^ { T } \right ) \, \rho \, d x \right ) W _ { + } \right ) & = N \lambda _ { \max } \left ( - W _ { + } ^ { T } C W _ { + } \right ) \\ & = N ( - \lambda _ { k } ( C ) ) & ( 3. 1 6 ) \\ & = N \lambda _ { m - k + 1 } ( - C ) \\ & = N \lambda _ { k ^ { \prime } } ( - C ), \\ \text{as required by Theorem 3.2. Set}$$

as required by Theorem 3.2. Set

$$\sigma _ { k ^ { \prime } } ^ { 2 } \, = \, \lambda _ { \max } \left ( \sum _ { j = 1 } ^ { N } W _ { + } ^ { T } ( L ^ { 2 } C ) W _ { + } \right ) \, = \, N L ^ { 2 } \lambda _ { \max } \left ( W _ { + } ^ { T } C W _ { + } \right ) \, = \, N L ^ { 2 } \lambda _ { 1 }. \ \ ( 3. 1 7 ) \\ \text{Theorem $3.2$ states}$$

Theorem 3.2 states

$$\mathbb { P } \left \{ \lambda _ { k ^ { \prime } } \left ( \sum _ { j = 1 } ^ { N } \left ( - \nabla _ { \mathbf x } f _ { j } \nabla _ { \mathbf x } f _ { j } ^ { T } \right ) \right ) \geq N \lambda _ { k ^ { \prime } } ( - C ) + N t \right \} \leq k \exp \left ( \frac { - ( N t ) ^ { 2 } } { 4 \sigma _ { k ^ { \prime } } ^ { 2 } } \right ) \quad ( 3. 1 8 ) \\ \text{when $N t \leq \sigma _ { k ^ { \prime } }^{2}/L^{2}$. Plug in the computed quantities with $t = - \varepsilon \lambda_{k ^ { \prime } } ( - C ) = \varepsilon \lambda_{k} ( C )$}$$

  when Nt ≤ σ 2 k ′ /L 2 . Plug in the computed quantities with t = -ελ k ′ ( -C ) = ελ k ( C ) to achieve the lower estimate (3.4). Note that the condition ε ≤ 1 ≤ λ /λ 1 k ′ allows us to apply Theorem 3.2.

Next we use this result to derive a lower bound on the number of gradient samples needed for relative accuracy of ε . Recall the definition of big omega notation that a = Ω( ) means b a ≥ cb for some positive constant c .

Corollary 3.3. Let κ k = λ /λ 1 k . Then for ε ∈ (0 , 1] ,

$$N \, = \, \Omega \left ( \frac { L ^ { 2 } \kappa _ { k } ^ { 2 } } { \lambda _ { 1 } \varepsilon ^ { 2 } } \log ( m ) \right )$$

implies | ˆ λ k -λ k | ≤ ελ k with high probability.

Proof . Starting with the upper estimate from Theorem 3.1, if

$$N \ \geq \ \frac { 4 L ^ { 2 } } { \lambda _ { k } \varepsilon ^ { 2 } } ( \beta + 1 ) \log ( m ) \ \geq \ \frac { 4 L ^ { 2 } } { \lambda _ { k } \varepsilon ^ { 2 } } ( \beta \log ( m ) + \log ( m - k + 1 ) ), \quad \ \ ( 3. 2 0 )$$

then

$$\mathbb { P } \left \{ \hat { \lambda } _ { k } \geq ( 1 + \varepsilon ) \lambda _ { k } \right \} \leq m ^ { - \beta }. \quad \quad \quad ( 3. 2 1 ) \\ \intertext { w e t i m a t e } \text{from Theorem $3 1$ if }$$

Similarly for the lower estimate from Theorem 3.1, if

$$N \, \geq \, \frac { 4 L ^ { 2 } \lambda _ { 1 } } { \lambda _ { k } ^ { 2 } \varepsilon ^ { 2 } } ( \beta + 1 ) \log ( m ) \, \geq \, \frac { 4 L ^ { 2 } \lambda _ { 1 } } { \lambda _ { k } ^ { 2 } \varepsilon ^ { 2 } } ( \beta \log ( m ) + \log ( k ) ), \quad \ ( 3. 2 2 )$$

then

$$\mathbb { P } \left \{ \hat { \lambda } _ { k } \leq ( 1 - \varepsilon ) \lambda _ { k } \right \} \leq m ^ { - \beta }. \quad \quad \quad ( 3. 2 3 ) \\. \text{ and taking} \sigma$$

Setting κ k = λ /λ 1 k and taking

$$N \, \geq \, ( \beta + 1 ) \frac { 4 L ^ { 2 } \kappa _ { k } ^ { 2 } } { \lambda _ { 1 } \varepsilon ^ { 2 } } \log ( m )$$

satisfies both conditions.

We can combine results from Golub and Van Loan [15, Chapter 8] with results from Tropp [30] to obtain an estimate of the distance between the subspace defined by the eigenvectors W 1 and the subspace defined by the eigenvectors ˆ W 1 . This requires a different matrix Bernstein inequality in the form of Theorem 6.1 from Tropp [30], which we restate below. When we apply the theorem, X j = ∇ x f j ∇ x f T j -C , that is, the random matrix samples are the deviation of the j th sampled gradient outer product from the matrix C .

Theorem 3.4 ( Matrix Bernstein: bounded case, Theorem 6.1 [30] ). Consider a finite sequence { X j } of independent, random, self-adjoint matrices with dimension n . Assume that

$$\mathbb { E } \left [ X _ { j } \right ] = 0 \quad a n d \quad \lambda _ { m a x } ( X _ { j } ) \leq R \quad a l m o s t \ s u r e l y.$$

Compute the norm of the total variance,

$$\sigma ^ { 2 } \coloneqq \left \| \sum _ { j } \mathbb { E } \left [ X _ { j } ^ { 2 } \right ] \right \|.$$

Then the following inequality holds for all τ ≥ 0 :

$$\cdot \quad \kappa \quad \kappa \quad \kappa \quad \kappa \quad \kappa \quad - \\ \mathbb { P } \left \{ \lambda _ { m a x } \left ( \sum _ { j } X _ { j } \right ) \geq \tau \right \} \leq \begin{cases} n \exp ( - 3 \tau ^ { 2 } / ( 8 \sigma ^ { 2 } ) ), & \tau \leq \sigma ^ { 2 } / R, \\ n \exp ( - 3 \tau / ( 8 R ) ), & \tau > \sigma ^ { 2 } / R. \end{cases} \\ \text{THEOREM} \ 3. 5. \ L e t \, \varepsilon > 0. \ A s s u m e \, \| \nabla _ { \mathbf x } f \| \leq L \ \text{for all} \, \mathbf x \in \mathcal { X }. \ \text{Define the variance} \, }$$

Theorem 3.5. Let ε &gt; 0 . Assume ‖∇ x f ‖ ≤ L for all x ∈ X . Define the variance and assume ν 2 &gt; 0 . Then

$$\nu ^ { 2 } \, = \, & \left \| \int \left ( \nabla _ { \mathbf x } f \, \nabla _ { \mathbf x } f ^ { T } - C \right ) ^ { 2 } \, \rho \, d x \right \|, \\ \Gamma ^ { 2 } > 0. \ T h e n & \quad \Gamma ^ { 2 } \, \cos _ { 2 } \sum _ { 2 } \sum _ { 2 } \sum _ { 2 } \sum _ { 2 } \sum _ { 2 } \sum _ { 2 } \sum _ { 2 } \sum _ { 2 } \sum _ { 2 } \sum _ { 2 } \sum _ { 2 } \sum _ { 2 } \sum _ { 2 } \sum _ { 2 } \sum _ { 2 } \sum _ { 2 } \sum _ { 2 } \sum _ { 2 } \sum _ { 2 } \sum _ { 2 } \sum _ { 2 } \sum _ { 2 } \sum _ { 2 } \sum _ { 2 } \sum _ { 2 } \sum _ { 2 } \sum _ { 2 } \sum _ { 2 } \sum _ { 2 } \sum _ { 2 } \sum _ { 2 } \sum _ { 2 } \sum _ { 2 } \sum _{ 2 } \sum _ { 2 } \sum _ { 2 } \sum _ { 2 } \sum _ { 2 } \sum _ { 2 } \sum _ { 2 } \sum _ { 2 } \sum _ { 2 } \sum _ { 2 } \sum _ { 2 } \sum _ { 2 } \sum _ { 2 } \sum _ { 2 } \sum _ { 2 } \sum _ { 2 } \sum _ { 2 } \sum _ { 2 } \sum _ { 2 } \sum _ { 2 } \sum _ { 2 } \sum _ { 2 } \sum _ { 2 } \sum _ { 2 } \sum _ { 2 } \sum _ { 2 } \sum _ { 2 } \sum _ { 2 } \sum _ { 2 } \sum _ { 2 } \sum _ { 2 } \sum _ { 2 } \sum _ { 2 } \sum _ { 2} \sum _ { 2 } \sum _ { 2 } \sum _ { 2 } \sum _ { 2 } \sum _ { 2 } \sum _ { 2 } \sum _ { 2 } \sum _ { 2 } \sum _ { 2 } \sum _ { 2 } \sum _ { 2 } \sum _ { 2 } \sum _ { 2 } \sum _ { 2 } \sum _ { 2 } \sum _ { 2 } \sum _ { 2 } \sum _ { 2 } \sum _ { 2 } \sum _ { 2 } \sum _ { 2 } \sum _ { 2 } \sum _ { 2 } \sum _ { 2 } \sum _ { 2 } \sum _ { 2 } \sum _ { 2 } \sum _ { 2 } \sum _ { 2 } \sum _ { 2 } \sum _ { 2 } \sum _ { 2 } \sum _ { 2 } \tau$$

$$& \mathbb { P } \left \{ \| \hat { C } - C \| \geq \varepsilon \| C \| \right \} \, \leq \, \begin{array} { c c } 2 m \exp \left ( \frac { - 3 N \lambda _ { 1 } ^ { 2 } \varepsilon ^ { 2 } } { 8 \nu ^ { 2 } } \right ), & \text{if $\varepsilon\leq\nu^{2}/(\lambda_{1}L^{2}$)}, \\ 2 m \exp \left ( \frac { - 3 N \lambda _ { 1 } \varepsilon } { 8 L ^ { 2 } } \right ), & \text{if $\varepsilon>\nu^{2}/(\lambda_{1}L^{2}$)}. \end{array} \\ & P r o o f \, \text{Observe that}$$

Proof . Observe that

Note that both

$$P r o o f. \text{ Observe that} \\ \mathbb { P } \left \{ \| \hat { C } - C \| \geq t \right \} = \mathbb { P } \left \{ \lambda _ { \max } ( \hat { C } - C ) \geq t \text{ or } \lambda _ { \max } ( C - \hat { C } ) \geq t \right \} \\ \leq \mathbb { P } \left \{ \lambda _ { \max } ( \hat { C } - C ) \geq t \right \} + \mathbb { P } \left \{ \lambda _ { \max } ( C - \hat { C } ) \geq t \right \} \\ = \mathbb { P } \left \{ \lambda _ { \max } \left ( \sum _ { j = 1 } ^ { N } \left ( \nabla _ { x } f _ { j } \nabla _ { x } f _ { j } ^ { T } - C \right ) \right ) \geq N t \right \} \\ + \mathbb { P } \left \{ \lambda _ { \max } \left ( \sum _ { j = 1 } ^ { N } \left ( C - \nabla _ { x } f _ { j } \nabla _ { x } f _ { j } ^ { T } \right ) \right ) \geq N t \right \}. \\ \text{Note that both} \\ \int \left ( \nabla _ { x } f \nabla _ { x } f ^ { T } - C \right ) \rho \, d x \, = \, \int \left ( C - \nabla _ { x } f \nabla _ { x } f ^ { T } \right ) \rho \, d x \, = \, 0.$$

$$\int \left ( \nabla _ { x } f \nabla _ { x } f ^ { T } - C \right ) \, \rho \, d x \, = \, \int \left ( C - \nabla _ { x } f \nabla _ { x } f ^ { T } \right ) \, \rho \, d x \, = \, 0. \quad \left ( 3. 2 8 \right ) \\ \text{Since $C$ is positive semidefinite and $\|\nabla_{x} f\| \leq L$,}$$

Since C is positive semidefinite and ‖∇ x f ‖ ≤ L ,

$$\lambda _ { \max } ( \nabla _ { \mathbf x } f \nabla _ { \mathbf x } f ^ { T } - C ) & = \max _ { \| \mathbf v \| = 1 } \mathbf v ^ { T } \left ( \nabla _ { \mathbf x } f \nabla _ { \mathbf x } f ^ { T } - C \right ) \mathbf v \\ & \leq \max _ { \| \mathbf v \| = 1 } \mathbf v ^ { T } \left ( \nabla _ { \mathbf x } f \nabla _ { \mathbf x } f ^ { T } \right ) \mathbf v \, \leq \, L ^ { 2 }. \\ \text{his bound also holds for } \lambda _ { \max } ( C - \nabla _ { \mathbf x } f \nabla _ { \mathbf x } f ^ { T } ), \, \since$$

This bound also holds for λ max ( C -∇ x f ∇ x f T ), since

$$\lambda _ { \max } ( C - \nabla _ { x } f \nabla _ { x } f ^ { T } ) & = \max _ { \| \mathbf v \| = 1 } \, \mathbf v ^ { T } \left ( C - \nabla _ { x } f \nabla _ { x } f ^ { T } \right ) \mathbf v \\ & \leq \max _ { \| \mathbf v \| = 1 } \, \mathbf v ^ { T } \, C \mathbf v \, \leq \, \| C \| \, \leq \, L ^ { 2 }.$$

Thus, the upper-bound R in Theorem 3.4 is L 2 . The variance parameter σ 2 is

$$\sigma ^ { 2 } = \left \| \left ( \sum _ { j = 1 } ^ { N } \int ( \nabla _ { \mathbf x } f _ { j } \nabla _ { \mathbf x } f _ { j } ^ { T } - C ) ^ { 2 } \, \rho \, d \mathbf x \right ) \right \| = N \left \| \int ( \nabla _ { \mathbf x } f \nabla _ { \mathbf x } f ^ { T } - C ) ^ { 2 } \, \rho \, d \mathbf x \right \| = N \, \nu ^ { 2 }.$$

Assume ε ≤ ν 2 / λ L ( 1 2 ). Then Nλ ε 1 ≤ Nν /L 2 2 , and we can apply the upper branch of Theorem 3.4 to the two terms at the end of (3.27) with τ = λ ε 1 = ‖ C ‖ ε , which produces the bound's upper branch in (3.26). Now assume ε &gt; ν 2 / λ L ( 1 2 ). Similarly, Nλ ε 1 ≥ Nν /L 2 2 , and we can apply the lower branch of Theorem 3.4 with the same t = λ ε 1 to (3.27) to produce the lower branch of (3.26).

This result leads to a lower bound on the number of samples needed for a small relative error in ˆ in the matrix 2-norm; compare the following corollary to Corollary C 3.3.

Corollary 3.6. Let ε &gt; 0 , and define

$$\delta \, = \, \max \left ( \frac { \nu ^ { 2 } } { \lambda _ { 1 } \varepsilon }, L ^ { 2 } \right ).$$

Then

$$N \, = \, \Omega \left ( \frac { \delta } { \lambda _ { 1 } \varepsilon } \, \log ( 2 m ) \right )$$

implies that ‖ ˆ C -C ‖ ≤ ε ‖ C ‖ with high probability.

Proof . We consider the two cases of ε from Theorem 3.5. Assume ε ≤ ν 2 / λ L ( 1 2 ), so that Nλ ε 1 ≤ Nν /L 2 2 , we follow the reasoning in the proof of Corollary 3.3 with the upper branch of the bound in Theorem 3.5 to get

$$N \, \geq \, \frac { 8 } { 3 } ( \beta + 1 ) \frac { \nu ^ { 2 } } { \lambda _ { 1 } ^ { 2 } \varepsilon ^ { 2 } } \log ( 2 m ).$$

Similarly, if ε &gt; ν 2 / λ L ( 1 2 ), then the lower branch from the bound in Theorem 3.5 produces

$$N \, \geq \, \frac { 8 } { 3 } ( \beta + 1 ) \frac { L ^ { 2 } } { \lambda _ { 1 } \varepsilon } \log ( 2 m ).$$

Using δ from (3.32) chooses the larger lower bound between (3.34) and (3.35).

We can combine Corollary 3.6 with [15, Corollary 8.1.11] to control the error in the estimated subspace defined by ˆ W 1 . We quantify this error by the distance between the subspace defined by W 1 and the subspace defined by ˆ W 1 . Recall the definition of the distance between subspaces [28],

$$\text{dist} ( \text{ran} ( W _ { 1 } ), \text{ran} ( \hat { W } _ { 1 } ) ) \, = \, \| W _ { 1 } W _ { 1 } ^ { T } - \hat { W } _ { 1 } \hat { W } _ { 1 } ^ { T } \| \, = \, \| W _ { 1 } ^ { T } \hat { W } _ { 2 } \|. \quad ( 3. 3 6 )$$

Corollary 3.7. Let ε &gt; 0 be such that

$$\varepsilon \, \leq \, \min \left ( 1, \left ( \lambda _ { n } - \lambda _ { n + 1 } \right ) / ( 5 \lambda _ { 1 } ) \right ),$$

and choose N according to Corollary 3.6. Then

$$\text{ dist} ( \tilde { W } _ { 1 } ), \, \tilde { \tilde { N } } _ { 1 } ) ) \, \leq \, \frac { 4 \lambda _ { 1 } \varepsilon } { \lambda _ { n } - \lambda _ { n + 1 } }, \, \text{ \quad \ \ } ( 3. 3 8 ) \, \text{ \quad \ } \text{ \quad \ } \text{ \quad \ } \text{ \quad \ } \text{ \quad \ } \text{ \quad \ } \text{ \quad \ } \text{ \quad \ } \text{ \quad \ } \text{ \quad \ } \text{ \quad \ } \text{ \quad \ } \text{ \quad \ } \text{ \quad \ } \text{ \quad \ } \text{ \quad \ } \text{ \quad \ } \text{ \quad \ } \text{ \quad \ } \text{ \quad \ } \text{ \quad \ } \text{ \quad \ } \text{ \quad \ } \text{ \quad \ } \text{ \quad \ } \text{ \quad \ } \text{ \quad \ } \text{ \quad \ } \text{ \quad \ } \text{\quad \ } \text{ \quad \ } \text{ \quad \ } \text{ \quad \ } \text{ \quad \ } \text{ \quad \ } \text{ \quad \ } \text{ \quad \ } \text{ \quad \ } \text{ \quad \ } \text{ \quad \ } \text{ \quad \ } \text{ \quad \ } \text{ \quad \ } \text{ \quad \ } \text{ \quad \ } \text{ \quad \ } \text{ \quad \ } \text{ \quad \ } \text{ \quad \ } \text{ \quad \ } \text{ \quad \ } \text{ \quad \ } \text{ \quad \ } \text{ \quad \ } \text{ \quad \ } \text{ \quad \ } \text{ \quad \ } \text{ \quad \ }$$

with high probability.

Proof . Let E = ˆ C -C . For ε in (3.37) with N chosen according to Corollary 3.6, we have

$$\| E \| \, \leq \, \varepsilon \| C \| \, = \, \varepsilon \lambda _ { 1 } \, \leq \, ( \lambda _ { n } - \lambda _ { n + 1 } ) / 5, \, \text{ \quad \ \ } ( 3. 3 9 )$$

with high probability. Under this condition on ‖ E ‖ , [15, Corollary 8.1.11] states

$$\text{dist} ( \text{ran} ( W _ { 1 } ), \text{ran} ( \hat { W } _ { 1 } ) ) \, \leq \, \frac { 4 \| E \| } { \lambda _ { n } - \lambda _ { n + 1 } } \, \leq \, \frac { 4 \lambda _ { 1 } \varepsilon } { \lambda _ { n } - \lambda _ { n + 1 } }, \quad \ \ ( 3. 4 0 )$$

as required.

Corollary 3.7 shows that control of the eigenvalues implies control of the subspace generated by the eigenvectors. However, the error in the estimated subspace is inversely proportional to the gap between the smallest eigenvalue associated with the active subspace and the largest eigenvalue associated with the inactive subspace. This implies, for example, if the gap between the second and third eigenvalues is larger than the gap between the first and second, then estimates of a two-dimensional active subspace are more accurate than estimates of a one-dimensional active subspace.

3.1. Approximate gradients. Many modern simulations have subroutines for estimating gradients with, e.g., adjoint methods [4, 3] or algorithmic differentiation [16]. However, legacy codes or simulations that couple multiple codes might not have such gradient capabilities. When there is no subroutine for gradients, finite difference approximations may suffice when m is not too large and f is neither too expensive nor too noisy. Recent work characterizes the gradient when function evaluations contain noise [22, 23].

Next, we extend the bounds on errors in the estimated eigenpairs to the case when the gradients are computed with some error. The gradient error model we analyze depends on a parameter that controls the amount of error. Let g x ( ) denote the approximate gradient computed at x ∈ X . We assume that

$$\lambda _ { \lambda } \quad \cdot \quad \cdot \\ \lambda _ { \lambda } \quad \| g ( x ) - \nabla _ { x } f ( x ) \| \, \leq \, \sqrt { m } \gamma _ { h }, \ \ x \in \mathcal { X },$$

where

$$\lim _ { h \to 0 } \gamma _ { h } \, = \, 0.$$

The parameter h may be a finite difference parameter, the grid spacing in a continuous adjoint computation, or the solver tolerance for a discrete adjoint computation.

Define the symmetric positive semidefinite matrix G and its eigenvalue decomposition

$$G \, = \, \int _ { \ } g g ^ { T } \, \rho \, d x \, = \, U \Theta U ^ { T }, \ \Theta = \text{diag} ( \theta _ { 1 }, \dots, \theta _ { m } ), \quad \ \ ( 3. 4 3 ) \\, \,, \,, \,, \,. \quad. \quad, \quad. \quad. \quad.$$

and define its random sample approximation

$$\hat { G } \, = \, \frac { 1 } { N } \sum _ { j = 1 } ^ { N } g _ { j } g _ { j } ^ { T } \, = \, \hat { U } \hat { \Theta } \hat { U } ^ { T }, \quad \hat { \Theta } = \text{diag} ( \hat { \theta } _ { 1 }, \dots, \hat { \theta } _ { m } ), \quad \quad ( 3. 4 4 ) \\. \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \ \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \ \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \,$$

where g j = g x ( j ) for x j drawn independently from ρ . With these quantities defined, we have the following lemma.

Lemma 3.8. Let ‖∇ x f ‖ ≤ L for all x ∈ X . The norm of the difference between ˆ C and ˆ G is bounded by

$$\hat { \ \ } \| \hat { C } - \hat { G } \| \, \leq \, ( \sqrt { m } \gamma _ { h } + 2 L ) \sqrt { m } \gamma _ { h }.$$

Proof . Let g = g x ( ) and ∇ x f = ∇ x f ( x ). First observe

$$\cdot \quad \sim \quad \sim \quad \sim \quad \sim \quad \sim \quad \sim \quad \sim \quad \sim \\ \| g + \nabla _ { x } f \| & = \| g - \nabla _ { x } f + 2 \nabla _ { x } f \| \leq \| g - \nabla _ { x } f \| + 2 \| \nabla _ { x } f \| \leq \sqrt { m } \gamma _ { h } + 2 L. \quad ( 3. 4 6 ) \\ \cdot \quad \quad \quad$$

Next,

$$\| g g ^ { T } - \nabla _ { \mathbf x } f \nabla _ { \mathbf x } f ^ { T } \| & = \frac { 1 } { 2 } \| ( g + \nabla _ { \mathbf x } f ) ( g - \nabla _ { \mathbf x } f ) ^ { T } + ( g - \nabla _ { \mathbf x } f ) ( g + \nabla _ { \mathbf x } f ) ^ { T } \| \\ & \leq \| ( g + \nabla _ { \mathbf x } f ) ( g -$$

Then,

$$\| \hat { G } - \hat { C } \| & = \left \| \frac { 1 } { N } \sum _ { j = 1 } ^ { N } \mathfrak { g } _ { j } ^ { T } - \frac { 1 } { N } \sum _ { j = 1 } ^ { N } \nabla _ { \mathbf x } f _ { j } \nabla _ { \mathbf x } f _ { j } ^ { T } \right \| \\ & \leq \frac { 1 } { N } \sum _ { j = 1 } ^ { N } \| \mathfrak { g } _$$

We combine Lemma 3.8 with Corollary 3.3 to study the error in the eigenvalues of the random sample estimate with approximate gradients.

Theorem 3.9. For ε ∈ (0 , 1] , if N is chosen as (3.19) , then the difference between λ k in (2.2) and the eigenvalue ˆ θ k from (3.44) is bounded as

$$\mu _ { \mu } \mu _ { \mu } \mu _ { \mu } \mu _ { \mu } \mu _ { \mu } \mu _ { \mu } \mu _ { \mu } \mu _ { \mu } \mu _ { \mu } \mu _ { \mu } \mu _ { \mu } \mu _ { \mu } \mu _ { \mu } \mu _ { \mu } \mu _ { \mu } \mu _ { \mu } \mu _ { \mu } \mu _ { \mu } \mu _ { \mu } \mu _ { \mu } \mu _ { \mu } \mu _ { \$$

with high probability.

Proof . Observe that

$$| \lambda _ { k } - \hat { \theta } _ { k } | \leq | \lambda _ { k } - \hat { \lambda } _ { k } | + | \hat { \lambda } _ { k } - \hat { \theta } _ { k } |. \quad \quad \quad ( 3. 5 0 )$$

Apply Corollary 3.3 to the first term. The second term follows from [15, Corollary 8.1.6] combined with Lemma 3.8, since

$$& | \hat { \theta } _ { k } - \hat { \lambda } _ { k } | = | \lambda _ { k } ( \hat { G } ) - \lambda _ { k } ( \hat { C } ) | \leq \| \hat { G } - \hat { C } \| \leq \sqrt { m } \gamma _ { h } ( \sqrt { m } \gamma _ { h } + 2 L ). \quad \ \ ( 3. 5 1 ) \\ - &$$

The bias in the finite sample eigenvalue estimates using approximate gradients goes to zero at the same rate as the error in the approximate gradient. Next we attend to the error in the active subspace computed with Monte Carlo and approximate gradients.

Theorem 3.10. Choose ε &gt; 0 such that

$$\varepsilon \, < \, \min \left ( 1, \, \frac { \lambda _ { n } - \lambda _ { n + 1 } } { 5 \lambda _ { 1 } }, \, \frac { \lambda _ { n } - \lambda _ { n + 1 } } { \lambda _ { n } + \lambda _ { n + 1 } } \right ) \, \quad \quad ( 3. 5 2 ) \\. \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdots$$

Choose N so that it satisfies both (3.33) and (3.19) with k = n +1 . Choose h small enough so that

$$\sqrt { m } \gamma _ { h } ( \sqrt { m } \gamma _ { h } + 2 L ) \leq \frac { ( 1 - \varepsilon ) \lambda _ { n } - ( 1 + \varepsilon ) \lambda _ { k + 1 } } { 5 }. \quad \quad ( 3. 5 3 )$$

Then

$$\ s u n d \quad & \quad \text{dist} ( \text{ran} ( \hat { U } _ { 1 } ), \, \text{ran} ( W _ { 1 } ) ) \, \leq \, \frac { 4 \sqrt { m } \gamma _ { h } ( \sqrt { m } \gamma _ { h } + 2 L ) } { ( 1 - \varepsilon ) \lambda _ { n } - ( 1 + \varepsilon ) \lambda _ { n + 1 } } \, + \, \frac { 4 \lambda _ { 1 } } { \lambda _ { n } - \lambda _ { n + 1 } }, \quad ( 3. 5 4 ) \\ \dots \, \text{and} \, \, \iota \colon \, \iota \dots \, \dots \, \iota \colon \iota \dots.$$

with high probability.

Proof . The conditions on N and ε imply | ˆ λ n +1 -λ n +1 | ≤ ελ n +1 with high probability due to Corollary 3.3. Examining (3.19), we see that if N is large enough to estimate λ n +1 , then N is large enough to estimate λ n , so | ˆ λ n -λ n | ≤ ελ n with high probability, too. Then

$$\lambda _ { n } - \lambda _ { n + 1 } & = | \lambda _ { n } - \lambda _ { n + 1 } | \\ & \leq | \lambda _ { n } - \hat { \lambda } _ { n } | + | \hat { \lambda } _ { n + 1 } - \lambda _ { n + 1 } | + ( \hat { \lambda } _ { n } - \hat { \lambda } _ { n + 1 } ) \quad \quad ( 3. 5 5 ) \\ & \leq \varepsilon \lambda _ { n } + \varepsilon \lambda _ { n + 1 } + ( \hat { \lambda } _ { n } - \hat { \lambda } _ { n + 1 } ), \\ \ h i b h \ n r o b i l i t v \text{.\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\} }$$

with high probability. Rearranging this inequality yields

$$\hat { \lambda } _ { n } - \hat { \lambda } _ { n + 1 } \, \geq \, ( 1 - \varepsilon ) \lambda _ { n } - ( 1 + \varepsilon ) \lambda _ { n + 1 }.$$

This relates the gap between the eigenvalue estimates to the gap between the true eigenvalues. The condition on ε ensures that

$$( 1 - \varepsilon ) \lambda _ { n } - ( 1 + \varepsilon ) \lambda _ { n + 1 } \, > \, 0.$$

Next,

$$\text{ dist} ( \hat { U } _ { 1 } ), \, \text{ ran} ( \hat { W } _ { 1 } ) ) \, \leq \, \text{ dist} ( \text{ran} ( \hat { U } _ { 1 } ), \, \text{ ran} ( \hat { W } _ { 1 } ) ) + \text{dist} ( \text{ran} ( \hat { W } _ { 1 } ), \, \text{ ran} ( \hat { W } _ { 1 } ) ). \\ \text{ w } \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \, \sim \, \cdot \cdot \cdot \cdot \cdot \cdot \cdot \, \sim \, \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \,.$$

The second term on the right is bounded in Corollary 3.7 under the assumptions on N and ε . The condition (3.53) on h and (3.56) imply

$$\sqrt { m } \gamma _ { h } ( \sqrt { m } \gamma _ { h } + 2 L ) \leq \frac { \hat { \lambda } _ { n } - \hat { \lambda } _ { k + 1 } } { 5 }.$$

Then [15, Corollary 8.1.11] implies

$$\text{ dist} ( \hat { U } _ { 1 } ), \, \text{ran} ( \hat { W } _ { 1 } ) ) \, \leq \, \frac { 4 } { \hat { \lambda } _ { n } - \hat { \lambda } _ { n + 1 } } \| \hat { G } - \hat { C } \|. \quad \quad ( 3. 6 0 ) \\. \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdot \, \cdots$$

Combining this with (3.56) and the bound from Lemma 3.8 yields the result.

In summary, the eigenvalues and the active subspace approximated with Monte Carlo and approximate gradients are well-behaved. The error bounds include a term that goes to zero like the error in the approximate gradient and a term that behaves like the finite sample approximation with exact gradients. Note that the error bound on the subspace estimate depends on both the gap between λ n and λ n +1 and a smaller gap that depends on ε .

- 4. Practical approach to computation. The bounds we present in Section 3 provide a theoretical foundation for understanding the behavior of the Monte Carlo estimates. However, many of the quantities in the bounds may not be known a priori-such as the maximum norm of the gradient L and the true eigenvalues of the matrix C . In this section we offer a practical recipe guided by the insights from the

theory. We caution that the following approach, which relies on a nonparametric bootstrap, can perform poorly for badly-behaved functions. For example, one could be unlucky and sample the gradient in regions that are not representative of the gradient over the entire domain; the bootstrap uses only the N samples used to compute the eigenpair estimates. Also, errors in the gradients could produce poor approximations of the eigenvalues and subspaces; we show an example of this in Section 5.1. Nevertheless, we have used the following approach on several problems in practice to reveal low-dimensional structure in complex functions of several variables coming from engineering simulations [7, 8, 21].

The first objective is to estimate the eigenvalues and a measure of the estimates' variability from the finite samples. Suppose one wishes to estimate the first k eigenpairs from the matrix C . Practical considerations guide the choice of k . For example, if one wishes to build a response surface approximation of f on a low-dimensional domain, then five or six dimensions might be the most one can afford given the cost of computing f ( x ). Hence k might be seven or eight to allow the possibility of finding a gap that indicates a sufficiently low-dimensional approximation. If a gap is not present in the first k eigenvalues, then f may not be amenable to dimension reduction via active subspaces for the desired purpose.

We recommend choosing the number N of independent gradient samples as

$$N \, = \, \alpha \, k \, \log ( m ),$$

where α is a multiplier between 2 and 10. Taking at least k samples means that C is a sum of k rank-one matrices and thus has a rank of at most k . This allows the possibility of estimating k non-zero eigenvalues. The log( m ) term follows from the bounds in Theorem 3.1. The α between 2 and 10 is an ad hoc multiplicative factor that we have used on several problems. In principle, αk is meant to model the contribution from the unknown terms L κ , k , ν 2 , and λ 1 in (3.19) and (3.33). It is likely that the combination of these terms with the ε -1 is greater than 10 k for small ε . However, the Bernstein inequalities used to derive the lower bounds on N in Corollaries 3.3 and 3.6 are also conservative. One can also assess if N is large enough a posteriori by examining the bootstrap intervals described below.

We form ˆ C using the samples of the gradient as in (3.1), and then compute its eigenvalue decomposition. We expect that computing the full eigendecomposition is much cheaper than computing the gradient samples. A function of a thousand variables produces ˆ C with dimension thousand-by-thousand. Full eigendecompositions for matrices this size are computed in seconds on modern laptops.

We suggest computing bootstrap intervals for the eigenvalues, which involves computing the eigendecompositions of several matrices the size of ˆ . C The bootstrap creates replicates by (i) sampling with replacement from the set of gradient samples, (ii) computing the replicate ˆ C ∗ , and (iii) computing its eigenvalue decomposition. The collection of eigenvalue replicates is used to estimate bounds on the true eigenvalues. Efron and Tibshirani use the bootstrap to get empirical density functions of estimated eigenvalues from a covariance matrix in section 7.2 of their book [10]. Chapter 3 of Jolliffe's book [19] also comments on the bootstrap approach for estimating eigenvalues and eigenvectors of a covariance matrix from independent samples. The bootstrap estimates of the standard error and confidence intervals for the eigenvalues may be biased, but this bias decreases as the number N of samples increases. Since these estimates may be biased, we refer to them as bootstrap intervals instead of confidence intervals .

Corollary 3.7 says that the error in the estimated subspace depends inversely on the gap between the eigenvalues scaled by the largest eigenvalue. The key to accurately approximating the subspace is to look for gaps in the eigenvalues; this is consistent with standard perturbation theory for eigenvector computations [28]. For example, if there is a larger gap between the third and fourth eigenvalues than between the second and third, then estimates of the three-dimensional subspace are more accurate than estimates of the two-dimensional subspace. This contrasts with heuristics for choosing the dimension of the subspace in (i) model reduction based on the proper orthogonal decomposition [2] and (ii) dimension reduction based on principal component analysis [19]. In these cases, one chooses the dimension of the subspace by a threshold on the magnitude of the eigenvalues-e.g., so that the sum of retained eigenvalues exceeds some proportion of the sum of all eigenvalues. To accurately approximate the active subspace, the most important quantity is the spectral gap, which indicates a separation between scales. To tease out the spectral gap, plot the estimated eigenvalues and their respective upper and lower bootstrap intervals; a gap between subsequent intervals offers confidence of a spectral gap and, hence, the presence of an active subspace. In Section 5, we show several examples of such plots (Figures 5.1, 5.2, 5.4, 5.5).

One should also consider the intent of the dimension reduction when choosing the dimension of the active subspace. For example, if the goal is to approximate a function of m variables by a surrogate model of n variables-as in [7]-then one may be limited to n 's small enough to permit surrogate construction. Suppose the largest n one is willing to use is n = n max = 5, but there is no gap between consecutive eigenvalues from λ 1 to λ 6 . Then subspace-based dimension reduction may not be an appropriate tool, and one should consider searching for other types of exploitable structure in the model.

Assuming we have chosen n , we wish to study the variability in the active subspace due to finite sampling; we again turn to the bootstrap. In particular, for each replicate ˆ W ∗ of the eigenvectors, we compute dist(ran( ˆ W 1 ) , ran( ˆ W ∗ 1 )). One can examine the bootstrap intervals of this quantity to study the stability of the subspace. Recall that the distance between subspaces is bounded above by 1, so a bootstrap interval whose values are close to 1 indicates a poorly approximated active subspace. Figures 5.1, 5.3, 5.4, and 5.5 show examples of plotting this metric for the stability of the subspace; the first two Figures also compare the measure of stability to the true error in the active subspace.

- 4.1. A step-by-step procedure. We summarize the practical approach to approximating the active subspace with bootstrap intervals. What follows is a modification of the procedure outlined at the beginning of Section 3 including our suggestions for parameter values. This procedure assumes the user has decided on the number k of eigenvalues to examine.
- 1. Choose N = αk log( m ), where α is a multiplier between 2 and 10, and choose N boot between 100 and 10000.
- 2. Draw N samples { x j } independently from ρ . For each x j , compute ∇ x f j = ∇ x f ( x j ).
- 3. Compute

$$\hat { C } \, = \, \frac { 1 } { N } \sum _ { j = 1 } ^ { N } ( \nabla _ { \times } f _ { j } ) ( \nabla _ { \times } f _ { j } ) ^ { T } \, = \, \hat { W } \hat { \Lambda } \hat { W } ^ { T }.$$

- 4. Bootstrap: For i from 1 to N boot , let /lscript i 1 , . . . , /lscript i N be N integers drawn randomly from { 1 , . . . , N } with replacement, and compute

$$\hat { C } _ { i } ^ { * } \, = \, \frac { 1 } { N } \sum _ { j = 1 } ^ { N } ( \nabla _ { \mathfrak x } f _ { \ell _ { j } ^ { i } } ) ( \nabla _ { \mathfrak x } f _ { \ell _ { j } ^ { i } } ) ^ { T } \, = \, \hat { W } _ { i } ^ { * } \, \hat { \Lambda } _ { i } ^ { * } \, \left ( \hat { W } _ { i } ^ { * } \right ) ^ { T }. \quad \quad ( 4. 3 ) \\ \intertext { n........................................................................................................................................................................................................ }$$

The asterisk denotes a bootstrap replicate. Then compute the subspace distance

$$d _ { i } ^ { * } \, = \, \text{dist} ( \text{ran} ( \hat { W } ), \, \text{ran} ( \hat { W } _ { i } ^ { * } ) ).$$

## 5. Compute the intervals

$$\left [ \min _ { i } \hat { \lambda } _ { j, i } ^ { * }, \max _ { i } \hat { \lambda } _ { j, i } ^ { * } \right ], \quad j = 1, \dots, k \quad \quad ( 4. 5 ) \\ \intertext { i t h d i o m o l f o m } \i t h d i o m o l f o m \hat { \lambda } ^ { * } \ i n / A 2 \diagdown \, \Lambda \iota \i o m u t h o m o l f o m$$

- 6. Plot the eigenvalue bootstrap intervals and look for large gaps. Choose the dimension n of the active subspace corresponding to the largest eigenvalue gap. If there is no perceivable gap, then an active subspace may not be present in the first k -1 dimensions.

where ˆ λ ∗ j,i is the j th diagonal from Λ ˆ ∗ i in (4.3). Also compute the mean, minimum, and maximum from the set { d ∗ i } to estimate the subspace error.

ˆ is small enough

A few comments are in order. First, we assume the dimension m of C so that the eigendecompositions of ˆ C and its bootstrap replicates are much cheaper than the samples of the gradient. Such is the case when m is in the thousands (i.e., f depends on thousands of input variables), and f and ∇ x f come from an expensive engineering simulation. Second, we choose the bootstrap to examine the variability because we assume that sampling more gradients is not feasible. If this is not the case, i.e., if one can cheaply evaluate many more gradient samples, then one can compute Monte Carlo estimates and central limit theorem confidence intervals of the eigenvalues in place of the bootstrap estimates. Lastly, we note that the elements of C are multivariate integrals. If m is small enough (2 or 3) and evaluating ∇ x f is cheap enough, then more accurate numerical quadrature rules may perform better than the random sampling, i.e., greater accuracy for fewer samples. However, practical error estimates are more difficult to compute, since the error is due to bias instead of variance.

- 5. Experiments. We apply the procedures described in Section 4 to two models: (i) a quadratic function and (ii) a linear functional of the solution of a parameterized PDE. The quadratic model is simple enough to analytically derive the eigenpairs of the active subspace for thorough evaluation of the method. We study the same PDE model in [7, Section 5]. Gradients are available through adjoints, but the true active subspaces are not available. We support efforts for reproducible research [20, 9]; codes for the experiments in this section can be found at

https://www.cs.purdue.edu/homes/dgleich/codes/compute-asm/compute-asm-code.tar.gz The PDE example uses the Random Field Simulation code ( http://www.mathworks.com/matlabcentral/fileexchange/27

as well as the MATLAB PDE Toolbox.

- 5.1. A quadratic model. Consider a quadratic function of m = 10 variables,

$$f ( x ) \, = \, \frac { 1 } { 2 } x ^ { T } A x, \ \ x \in [ - 1, 1 ] ^ { 1 0 },$$

where A is symmetric and positive definite. We take ρ = 2 -10 on the hypercube [ -1 1] , 10 and zero elsewhere. The gradient is ∇ x f ( x ) = A x , so

$$C \, = \, A \left ( \int x x ^ { T } \, \rho \, d x \right ) A ^ { T } \, = \, \frac { 1 } { 3 } A ^ { 2 }.$$

The eigenvalues of C are the squared eigenvalues of A divided by 3, and the eigenvectors of C are the eigenvectors of A .

We study three different A 's constructed from three choices for the eigenvalues: (1) exponential decay with a constant rate, (2) like the first but with a larger gap between the first and second eigenvalue, and (3) like the first with a larger gap between third and fourth eigenvalue. The three cases of eigenvalues for A are shown in the top row of Figure 5.1. Each A has the same eigenvectors, which we generate as an orthogonal basis from a random 10 × 10 matrix.

To estimate the eigenvalues, we choose N as in (4.1) with the multiplier α = 2 and k = 6 eigenvalues of interest, which yields N = 28 evaluations of the gradient. The middle row of Figure 5.1 shows the bootstrap intervals for the first six eigenvalues along with the true eigenvalues of C . The small bootstrap intervals suggest confidence in the estimates. The gaps are apparent in the last two cases. The bottom row of Figure 5.1 shows bootstrap intervals on the distance between the true k -dimensional active subspace and the subspace estimated with the N samples; the true distance is indicated by the circles. Notice that subspaces corresponding to the larger eigenvalue gap are much better approximated than the others. For example, the three-dimensional subspace is better approximated than the one- and two-dimensional subspaces for the third case.

Next we repeat the study using finite difference approximations of the gradient with step size h = 10 -1 , 10 -3 , and 10 -5 . The first of these step sizes is larger than would normally be used for such a model. We chose this large value to study the interplay between inaccurate gradients and the finite sample approximations of the eigenpairs. Figure 5.2 shows the true eigenvalues, their estimates, and the bootstrap intervals for all three cases and all three values of h ; the horizontal lines show the value of h . Eigenvalues that are smaller than h are estimated less accurately than those larger than h , which is not surprising since we are using first order finite differences. Also the gaps are much less noticeable in the estimates when finite difference parameter is not small enough to resolve the smaller eigenvalue in the pair defining the gap. In fact, this particular problem shows a large gap in the finite difference approximations when there is none in the true eigenvalues; see Figure 5.2(b,c) for examples of this phenomenon.

Figure 5.3 shows the distance between the true active subspace and the finite sample estimate with approximate gradients (circles). We use the bootstrap to estimate the error in the subspace as in Section 4. There is a strong bias in the estimates of the subspace error when the corresponding eigenvalues are not properly resolved. For instance, in Figure 5.3, the estimates of the error for subspaces of dimension 4 through 6 are biased for h = 10 -3 and significantly biased for h = 10 -1 . Compare this to the error in the last three eigenvalues for the smallest h = 10 -5 in Figure 5.2.

5.2. A parameterized PDE model. In previous work [7], we exploited the active subspace in the following parameterized PDE model to efficiently construct a kriging surface. Here we perform a more careful study of the variation in the active subspace estimated with finite samples of the gradient. Consider the following linear

Fig. 5.1: The top row shows the eigenvalues of the three choices for A . The second row shows the true and estimated eigenvalues along with the bootstrap intervals; eigenvalues are well approximated for all three cases. The third row shows the distance between the estimated subspace and the true subspace. In practice we do not have the true subspace, but we can estimate the distance with a bootstrap procedure as described in Section 4; the bootstrap intervals are shown, and the accuracy of the subspace estimates corresponds to the gaps in the eigenvalues of C .

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0545v2-figure-19393c86-47dd-48a7-bc7e-8f695455a515.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0545v2-figure-08740677-d289-445d-af62-96d0a129a383.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0545v2-figure-d2892e23-96d6-4cb3-8b20-9f6c3e0bb903.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0545v2-figure-096d6f5a-8ebc-412c-b716-2565f58f4e69.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0545v2-figure-3e50f340-89e6-4698-b7ea-25b1f6a77cfa.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0545v2-figure-d3ff13a8-67f4-465d-ae36-b4857850b588.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0545v2-figure-d35a3288-58ba-4c70-9ef0-903d9c30dd98.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0545v2-figure-b89da504-1b3e-4eae-95e1-2670275d7fab.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0545v2-figure-9d15abbd-2046-41ef-b158-07038c85081a.jpg)

elliptic PDE with parameterized, variable coefficients. Let u = u ( s x , ) satisfy

$$- \, \nabla _ { s } \cdot ( a \, \nabla _ { s } u ) \, = \, 1, \quad s \in [ 0, 1 ] ^ { 2 }.$$

We set homogeneous Dirichlet boundary conditions on the left, top, and bottom of the spatial domain [0 , 1] 2 ; denote this boundary by Γ . 1 The right side of the spatial domain, denoted Γ , has a homogeneous Neumann boundary condition. 2 The log of

Fig. 5.2: Eigenvalues, estimates, and bootstrap intervals using finite difference gradients with h = 10 -1 (top row), h = 10 -3 (middle row), and h = 10 -5 (bottom row). The horizontal black lines indicate the value of h in each plot. In general, estimates of eigenvalues smaller than h are less accurate than those larger than h .

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0545v2-figure-9f1ddf53-4dc6-4c94-9639-a2f2f9c6191f.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0545v2-figure-3aed872a-e9a1-47a0-adb4-3ecbcb701eb6.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0545v2-figure-2bc5c21c-d557-4ff6-982f-a865afb065f9.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0545v2-figure-92f94583-35a9-4224-8119-286ab11c6229.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0545v2-figure-e59e60dd-1c8e-4c91-a33a-3c04c7afeb25.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0545v2-figure-1b365fb8-ac17-4104-8bbf-8b2b4f14ead4.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0545v2-figure-80510847-1313-416c-ad01-f26e4f9c1af6.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0545v2-figure-7084bb0b-4317-4789-b586-bcdb848d102b.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0545v2-figure-a1e864cc-1c65-4023-bcc3-adc06a1cd4ad.jpg)

the coefficients a = a ( s x , ) is given by a truncated Karhunen-Loeve-type expansion

$$\log ( a ( s, x ) ) \, = \, \sum _ { i = 1 } ^ { m } x _ { i } \, \gamma _ { i } \, \phi _ { i } ( s ), & & ( 5. 4 ) \\ \text{indonond} \, \text{dontroll} \, \text{distributed standard normal random vari}$$

where the x i are independent, identically distributed standard normal random variables, and the { φ i ( s ) , γ i } are the eigenpairs of the correlation operator

$$\mathcal { C } ( s, t ) \, = \, \exp \left ( - \beta ^ { - 1 } \left \| s - t \right \| _ { 1 } \right ). \quad \quad \quad ( 5. 5 ) \\ \text{ality of the active subspace approximation for two correlation lengths, } \cdots \, \cdots \, \cdots \, \cdots \,.$$

We study the quality of the active subspace approximation for two correlation lengths, β = 1 and β = 0 01. These correspond to . long and short correlation lengths, respec-

Fig. 5.3: The distance between the true active subspace and its finite sample approximation along with bootstrap intervals for h = 10 -1 (top row), h = 10 -3 (middle row), and h = 10 -5 (bottom row). The subspaces are very poorly approximated when the finite difference step size is not small enough to resolve the eigenvalues corresponding the subspaces; compare to Figure 5.2. However, subspaces with a larger associated eigenvalue gap are generally approximated better than others.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0545v2-figure-be18f937-7f52-44fd-85a0-16fafdf0f982.jpg)

tively, for the random field defining the log of the coefficients. We choose a truncation of the field m = 100, which implies that the parameter space X = R 100 with ρ a standard Gaussian density function. Define the linear function of the solution

$$f ( \mathbf x ) \, = \, \frac { 1 } { | \Gamma _ { 2 } | } \int _ { \Gamma _ { 2 } } u ( \mathbf s, \mathbf x ) \, d \mathbf s. & & ( \mathbf 5. 6 ) \\ \mathcal { H } \, \cdots \, \mathcal { H } \, \cdots \, \mathcal { H } \, \cdots \, \mathcal { H } \, \cdots \, \mathcal { H } \, \cdots \, \mathcal { H } \, \cdots \, \mathcal { H } \, \cdots \, \mathcal { H } \, \cdots \, \mathcal { H } \, \cdots \, \mathcal { H } \, \cdots \, \mathcal { H } \, \cdots \, \mathcal { H } \, \cdots \, \mathcal { H } \, \cdots \, \mathcal { H } \, \cdots \, \mathcal { H } \, \cdots \, \mathcal { H } \, \cdots \, \mathcal { W } \, \cdots \, \mathcal { H } \, \cdots \, \mathcal { W } \, \cdots \, \mathcal { W } \, \cdots \, \mathcal { W } \, \cdots \, \mathcal { W } \, \cdots \, \mathcal { W } \, \cdots \, \mathcal { W } \, \cdots \, \mathcal { W } \, \cdots \, \mathcal { W } \, \cdots \, \mathcal { W } \, \cdots \, \mathcal { W } \, \cdots \, \mathcal { W } \, \cdots \, \mathcal { W } \, \cdots \, \mathcal { W } \, \cdots \, \mathcal { W } \, \cdots \, \mathcal { W } \, \cdots \, \mathcal { W } \, \mathcal { W } \, \cdots \, \mathcal { W } \, \cdots \, \mathcal { W } \, \cdots \, \mathcal { W } \, \cdots \, \mathcal { W } \, \cdots \, \mathcal { W } \, \cdots \, \mathcal { W } \, \cdots \, \mathcal { W } \, \cdots \, \mathcal { W } \, \cdots \, \mathcal { W } \, \cdots \, \mathcal { W } \, \cdots \, \mathcal { W } \, \cdots \, \mathcal { W } \, \cdots \, \mathcal { W } \, \cdots \, \mathcal { W } \, \cdots \, \mathbb { W } \, \cdots \, \mathcal { W } \, \cdots \, \mathcal { W } \, \cdots \, \mathcal { W } \, \cdots \, \mathcal { W } \, \cdots \, \mathcal { W } \, \cdots \, \mathcal { W } \, \cdots \, \mathcal { W } \, \cdots \, \mathcal { W } \, \cdots \, \mathcal { W } \, \cdots \, \mathcal { W } \, \cdots \, \mathcal { W } \, \cdots \, \mathcal { W } \, \cdots \, \mathcal { W } \, \cdots \, \mathcal { W } \, \cdots \, \mathcal { W } \cdots \, \mathcal { W } \, \cdots \, \mathcal { W } \, \cdots \, \mathcal { W } \, \cdots \, \mathcal { W } \, \cdots \, \mathcal { W } \, \cdots \, \mathcal { W } \, \cdots \, \mathcal { W } \, \cdots \, \mathcal { W } \, \cdots \, \mathcal { W } \, \cdots \, \mathcal { W } \, \cdots \, \mathcal { W } \, \cdots \, \mathcal { W } \, \cdots \, \mathcal { W } \, \cdots \, \mathcal { W } \, \cdots \, \mathcal { W } \, \cdots \$$

This is the quantity of interest from the model (more precisely, its approximation with a finite element method). Given a value for the input parameters x , we discretize the PDE with a standard linear finite element method using MATLAB's PDE Toolbox. The discretized domain has 34320 triangles and 17361 nodes; the eigenfunctions φ i

from (5.4) are approximated on this mesh. We compute the gradient of the quantity of interest (5.6) using a discrete adjoint formulation. Further details appear in our previous work [7].

The top row of Figure 5.4 shows the estimates of the eigenvalues of C along with the bootstrap intervals for β = 1 in (5.5). The gap between the first and second eigenvalues is apparent and supported by the gap in the corresponding bootstrap intervals. We exploit this gap in [7] to construct an accurate univariate kriging surface of the active variable. The bottom row of Figure 5.4 shows the variance in the estimated subspace as computed with the bootstrap including the bootstrap intervals. The left column of Figure 5.4 uses the multiplier α = 2 when choosing the number N of gradient samples; the right column uses α = 10. Notice the overall decrease in both the range of the bootstrap interval and the subspace error as we include more samples. Figure 5.5 shows the identical study with the short correlation length β = 0 01 from . (5.5).

6. Summary and conclusions. Consider a scalar-valued function of several variables. The average outer product of the gradient with itself is the central matrix in the development of active subspaces for dimension reduction. The dominant eigenvectors define the directions along which input perturbations change the output more, on average. We have analyzed a Monte Carlo method for approximating this matrix and its eigenpairs. We use recent theory developed for the eigenvalues of sums of random matrices to analyze the probability that the finite sample eigenvalue estimates deviate from the true eigenvalues, and we combine this analysis with results from matrix computations to derive results for the subspaces. We extend this analysis to quantities computed with samples of approximate gradients, e.g., finite differences. We also provide a practical computational approach that employs the bootstrap to reveal the error in the eigenvalues and the stability of the subspace.

Our analysis offers answers to the following important questions. First, how many gradient samples does one need for an accurate approximation of the first k eigenvalues? Precise theoretical bounds motivate a heuristic that chooses a number proportional to k times the log of the dimension m . Second, what can be said about the accuracy of the estimated subspace? The accuracy of the estimated subspace is directly related to gaps in the eigenvalues. Third, how does one judge the stability of the computed quantities? We propose to use bootstrap intervals for the eigenvalues and the stability of the subspace. Finally, how does this analysis change when gradients are not exact but approximate? Our theory shows that approximate gradients introduce a bias term in the error bounds that goes to zero as the approximate gradients become more accurate. The numerical examples suggest that this bias can produce inaccurate subspaces when the gradients are not well approximated.

Acknowledgments. The first author was partially supported by the U.S. Department of Energy Office of Science, Office of Advanced Scientific Computing Research, Applied Mathematics program under Award Number DE-SC-0011077. The second author would like to thank the Simon's Institute for Theory of Computing program on Big Data for the opportunity to learn about the randomized methods used in this paper and NSF CAREER award CCF-114975.

## REFERENCES

[1] Athanasios C Antoulas , Approximation of large-scale dynamical systems , vol. 6, Society for Industrial and Applied Mathematics, 2005.

Eigenvalues

Subspace Distance

-6

10

-7

10

-8

10

-9

10

-10

10

-11

10

-12

10

-13

1

1

10

0

10

-1

10

-2

10

2

(a)

β

3

4

Index

= 1,

α

= 2

2

3

4

5

Subspace Dimension

(c)

3

4

Index

(b)

β

= 1,

α

= 10

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0545v2-figure-a94a9354-9205-4a51-951e-a55e9a4b43f5.jpg)

β

= 1,

α

= 2

Fig. 5.4: The top row shows estimates of the eigenvalues of C along with the bootstrap intervals for the quantity of interest (5.6) from the parameterized PDE model with the long correlation length β = 1 from (5.5). The bottom row shows the estimates and bootstrap intervals on the distance between the estimated active subspace and the true active subspace. The left column is computed with the multiplier α = 2 when choosing N ; the right column uses α = 10. The gap between the first and second eigenvalue is significant as judged by the gap between the bootstrap intervals.

- [2] G Berkooz, P Holmes, and J L Lumley , The proper orthogonal decomposition in the analysis of turbulent flows , Annual Review of Fluid Mechanics, 25 (1993), pp. 539-575.
- [3] Alfio Borzi and Volker Schulz , Computational Optimization of Systems Governed by Partial Differential Equations , SIAM, 2012.
- [4] Arther E. Bryson and Yu-Chi Ho , Applied Optimal Control: Optimization, Estimation, and Control , Hemisphere Publishing Corportation, 1975.
- [5] Hans-Joachim Bungartz and Michael Griebel , Sparse grids , Acta Numerica, 13 (2004), pp. 147-269.
- [6] Russel E. Caflisch , Monte carlo and quasi-monte carlo methods , Acta Numerica, 7 (1998), pp. 1-49.
- [7] P. Constantine, E. Dow, and Q. Wang , Active subspace methods in theory and practice: Applications to kriging surfaces , SIAM Journal on Scientific Computing, 36 (2014),

5

Est

BI

BI

Est

6

6

Eigenvalues

-6

10

-7

10

-8

10

-9

10

-10

10

-11

10

-12

10

-13

10

1

2

5

Est

BI

6

Fig. 5.5: The top row shows estimates of the eigenvalues of C along with the bootstrap intervals for the quantity of interest (5.6) from the parameterized PDE model with the short correlation length β = 0 01 from (5.5). The bottom row shows the estimates . and bootstrap intervals on the distance between the estimated active subspace and the true active subspace. The left column is computed with the multiplier α = 2 when choosing N ; the right column uses α = 10. The gap between the first and second eigenvalue is significant as judged by the gap between the bootstrap intervals.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0545v2-figure-fafa1f9e-fee9-46d7-9b0f-42cfd02f4d00.jpg)

## pp. A1500-A1524.

- [8] Paul G Constantine, Brian Zaharatos, and Mark Campanelli , Discovering an active subspace in a single-diode solar cell model , arXiv preprint arXiv:1406.7607, (2014).
- [9] David L Donoho, Arian Maleki, Inam Ur Rahman, Morteza Shahram, and Victoria Stodden , Reproducible research in computational harmonic analysis , Computing in Science &amp; Engineering, 11 (2009), pp. 8-18.
- [10] Bradley Efron and Robert J Tibshirani , An Introduction to the Bootstrap , vol. 57, CRC press, 1994.
- [11] Massimo Fornasier, Karin Schnass, and Jan Vybiral , Learning functions of few arbitrary linear parameters in high dimensions , Foundations of Computational Mathematics, 12 (2012), pp. 229-262.
- [12] Kenji Fukumizu and Chenlei Leng , Gradient-based kernel dimension reduction for regres-

sion , Journal of the American Statistical Association, 109 (2014), pp. 359-370.

- [13] Alex Gittens and Michael W. Mahoney , Revisiting the Nystr¨m method for o improved large-scale machine learning. , in ICML, vol. 28, 2013, pp. 567-575.
- [14] Alex Gittens and Joel A Tropp , Tail bounds for all eigenvalues of a sum of random matrices , arXiv preprint arXiv:1104.4513, (2011).
- [15] Gene H Golub and Charles F Van Loan , Matrix Computations , Johns Hopkins University Press, 3rd ed., 1996.
- [16] Andreas Griewank , Evaluating Derivatives: Principles and Techniques of Algorithmic Differentiation , Society for Industrial and Applied Mathematics, 2000.
- [17] N. Halko, P. G. Martinsson, and J. A. Tropp , Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions , SIAM Rev., 53 (2011), pp. 217-288.
- [18] Marian Hristache, Anatoli Juditsky, Jorg Polzehl, and Vladimir Spokoiny , Structure adaptive approach for dimension reduction , The Annals of Statistics, 29 (2001), pp. 15371566.
- [19] Ian Jolliffe , Principal Component Analysis , Wiley Online Library, 2005.
- [20] Randall J LeVeque , Python tools for reproducible research on hyperbolic problems , Computing in Science &amp; Engineering, 11 (2009), pp. 19-27.
- [21] Trent W. Lukaczyk, Paul Constantine, Francisco Palacios, and Juan J. Alonso , Active Subspaces for Shape Optimization , American Institute of Aeronautics and Astronautics, 2014/02/16 2014.
- [22] Jorge J. Mor´ e and Stefan M. Wild , Estimating derivatives of noisy simulations , ACM Trans. Math. Softw., 38 (2012), pp. 19:1-19:21.
- [23] Jorge J. Mor´ e and Stefan M. Wild , Do you trust derivatives or differences? , Journal of Computational Physics, 273 (2014), pp. 268 - 277.
- [24] Art B. Owen , Monte Carlo theory, methods and examples , 2013. http://statweb.stanford.edu/ ~ owen/mc/ .
- [25] Trent M. Russi , Uncertainty Quantification with Experimental Data and Complex System Models , PhD thesis, UC Berkeley, 2010.
- [26] A. Saltelli, M. Ratto, T. Andres, F. Campolongo, J. Cariboni, D. Gatelli, M. Saisana, and S. Tarantola , Global Sensitivity Analysis: The Primer , John Wiley &amp; Sons, 2008.
- [27] A.M. Samarov , Exploring regression structure using nonparametric functional estimation , Journal of the American Statistical Association, 88 (1993), pp. 836-847.
- [28] G. Stewart , Error and perturbation bounds for subspaces associated with certain eigenvalue problems , SIAM Review, 15 (1973), pp. 727-764.
- [29] Alex Townsend and Lloyd N. Trefethen , Continuous analogues of matrix factorizations , Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, 471 (2014), p. 20140585.
- [30] Joel A Tropp , User-friendly tail bounds for sums of random matrices , Foundations of Computational Mathematics, 12 (2012), pp. 389-434.
- [31] Yingcun Xia , A constructive approach to the estimation of dimension reduction directions , The Annals of Statistics, 35 (2007), pp. 2654-2690.
Paul Constantine David Gleich
Computing active subspaces with Monte Carlo
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.