
\section{Introduction}
	
	Simulation-based computer experiments are pervasive in diverse fields
	ranging from engineering to social science. Some of these simulations are
	time--and resource--consuming to such an extent that carefully selecting
	the experiments to be run is a keystone concern. A standard technique in
	this context is to replace the expensive simulator by an inexpensive
	surrogate (e.g., a quadratic response surface, a radial basis function, or
	a Gaussian process (GP)). Such approaches have been shown to be efficient,
	with the downside that performance drops sharply when the number of
	variables determining the experiment becomes more than a few dozen.
	
	Many efforts deal with this ``curse of dimensionality" (as coined by
	\citep{Bellman2003}) from different points of view. \resp{1-3}{Global} sensitivity analysis
	is concerned with determining the relative importance of each variable
	with respect to a quantity of interest, using methods such as screening or
	measures of influence; we refer to the work of \citet{Iooss2015} for a
	review. Another option is to assume additional structure about the problem
	at hand, such as additivity. The idea is to decompose the multivariate
	function $f$ into a sum of functions each of fewer variables, decreasing
	the amount of data needed for reliable inference. For GPs, purely additive
	examples include those of \citet{Durrande2010} and \citet{Duvenaud2011},
	while combinations of variables are considered, for instance, by
	\citet{Durrande2013} and \citet{Sung2017}. Learning a latent space has
	also been attempted, for example with a generative topographic mapping
	\citep{Viswanath2011} or internally in the GP model \citep{Titsias2010}.
	Related to this latter method is detecting the presence of an
	\textit{active subspace} \citep{Constantine2015}, the framework in which
	this article operates (see Section~\ref{sec:asm} for a formal definition).
	The principle is to find the directions accounting for most of the
	variation in the quantity of interest and then to conduct a study with a
	surrogate on those main directions only. This approach has been shown to
	perform well on a variety of problems, is easy to interpret because it
	corresponds to learning important linear combinations of the inputs, and
	has a theoretical foundation. It additionally allows for some
	visualization when the number of directions appears (or is selected) to be
	small.
	
	The methodology typically consists of a random sampling of the
	high-dimensional design space to identify the active subspace, followed by
	a parameter study therein. While some guidelines for the primary stage
	exist, the methodology remains unsatisfying for several reasons. First, it
	requires gradient observations, which are not always available for legacy
	simulators or noisy problems \citep{LMW2019AN}. Relying on
	finite-difference approximations is an option--even for noisy simulators
	\citep{More2012}--but it has been shown to be less desirable
	\citep{Palar2018}. Directional derivatives may be used as well, from a
	low-rank matrix recovery perspective
	\citep{Djolonga2013,Constantine2015b}, but still demanding more costly
	observations if not directly available. As a result, replacing these
	derivative observations with those from a surrogate (e.g., a GP) is a
	standard approach \citep{Fukumizu2014,Othmer2016,Palar2018}. One must,
	however, still quantify the additional amount of uncertainty introduced.
	Second, unless a surrogate is used, the learned subspace cannot benefit
	from experiments run during the second stage because of the assumption of
	independent and identically distributed (iid) gradient sampling (because,
	for instance, optimization will typically involve function evaluations
	clustered in specific locations).
	
	In the GP literature, several authors have bypassed the two-stage approach
	by directly estimating the directions of influence, treating those as
	additional hyperparameters; see, for example, the work of
	\citet{Garnett2014} and \citet{Tripathy2016}.

	As a byproduct, this approach does not require gradient observations
	(although such observations can be accommodated) and does not require iid
	samples. \resp{}{Also, \citet{Marcy2017} proposed a sophisticated 
	fully Bayesian estimation framework of the dimension reduction
	matrix (including its dimension), which entails a challenging
	 optimization over the Grasmannian or Stiefold manifold.} 
 	Nevertheless, the proposed estimation procedures
	become much more complicated, requiring a variety of approximations. A
	major benefit of not requiring iid samples for active subspace estimation
	is to enable sequential learning procedures, that is, adding optimally
	informative points with respect to a quantity of interest.
	
	
	In this article, we bridge the gap between two-stage and direct estimation
	procedures for GPs, alleviating some of the limitations of recent active
	subspace methods. Our contributions can be summarized as follows:
	\begin{itemize}
		\item We prove that the classical active subspace estimator is a
		method of moments estimator for a certain class of GP models, formally
		establishing a link where similarities have previously been noted.
		\item We provide a closed-form expression for the matrix defining the
		active subspace for a GP for popular covariance kernels, rendering the
		Monte Carlo (MC) finite differencing on top of a surrogate
		unnecessary;
		\item Our closed-form solution allows sequential design to be carried
		out efficiently, allowing us to gain better estimates of the active
		subspace with fewer black-box evaluations, a feature especially important 
		when these evaluations are computationally expensive;
	
		\item We show how interval estimates may be calculated quickly via MC
		on the active subspace by propagating uncertainty from the GP
		hyperparameters.
	\end{itemize} 
	
	This paper is organized as follows. Section~\ref{sec:background} presents
	the active subspace approach as well as background on GPs.
	Section~\ref{sec:methodology} contains our methodological contributions,
	before illustrating them empirically in Section~\ref{sec:experiments}. We
	present our conclusions in Section~\ref{sec:discussion}.
	
	
	\section{Background}
	\label{sec:background}
	
	We first review the active subspace and GP paradigms.
	
	\subsection{Active Subspaces}
	
	\label{sec:asm} 
	
    The concept of \textit{active subspaces} \citep{Constantine2014} is
    appealing for fields from uncertainty quantification (UQ) to inverse
    problems. Informally, active subspace methods (ASMs) involve determining
    along which directions a scalar function of many variables changes a lot
    on average, and along which directions it is almost constant on average.
    The ideal family functions with which to describe the ASMs is that of
    \textit{ridge functions}, (not to be confused with ridge regression, an
    unrelated $\ell_2$ regularization technique), which are defined as
    functions that vary \textit{uniquely} in certain directions. To be
    precise, fix $\mathbf{x} \in \mathcal{X} \subseteq \mathbb{R}^m$ and $\mathbf{A} \in \mathbb{R}^{r
    \times m}$. Then a ridge function $f:\mathbb{R}^m \rightarrow \mathbb{R}$ is a
    function of the form
	\begin{equation}
	f(\mathbf{x}) = g(\mathbf{A} \mathbf{x})
	\label{ridgefunc}
	\end{equation}	
    for some function $g:\mathbb{R}^r \rightarrow \mathbb{R}$. Typically, $r < m$, so $g$
    acts on a lower dimension than $m$. By discovering $\mathbf{A}$, one
    exposes the low-dimensional structure in $f$. These functions underlie
    projection pursuit methods; see, for instance, the work of
    \citet{Friedman1981}.

    Ridge functions are constant along $\mathbf{A}$'s kernel. ASMs provide a
    framework for operating under looser assumptions: instead of the function
    being constant along $\mathbf{A}$'s kernel, ASMs stipulate that the
    directional derivatives in directions belonging to this subspace are
    significantly smaller than those belonging to $\mathbf{A}$'s range and
    make this assumption in expectation rather than absolutely.
	
	Rigorous definition of the active subspace requires us to consider the
	matrix $\mathbf{C}$, defined as the expected outer product of the gradient:
	\begin{equation}
	\mathbf{C} = \int_{\mathcal{X}} (\nabla f) (\nabla f)^\top d\mu.
	\label{Cdef}
	\end{equation}
    Here, $\mathcal{X}\subseteq \mathbb{R}^m$ is the domain (of interest) of $f$, and
    $\mu$ is an arbitrary measure. Little attention has been paid to
    specifying $\mu$ beyond the heuristic that the Lebesgue measure is in some
    sense the natural starting point for bounded $\mathcal{X}$, whereas the Gaussian
    is used for unbounded $\mathcal{X}$.
	
	In practice, $\mathbf{C}$ is estimated by using a simple Monte Carlo procedure:
	sample $\mathbf{x}_{i} \sim \mu$ for $i \in \{1, \ldots, M\}$, then calculate the
	gradient at each point $\mathbf{x}_i$ to obtain the estimator
	\begin{equation}
	\hat{\mathbf{C}} = \frac{1}{M}\sum_{i=1}^M (\nabla f(\mathbf{x}_i)) (\nabla f(\mathbf{x}_i))^\top.
	\label{Chatdef}
	\end{equation}
	The properties of this estimator, which we call the classical or MC
	estimator, have been analyzed by \citet{Constantine2015} and
	\citet{Holodnak2018}, along with guidelines for choosing $M$, the MC
	sample size.
	
    The matrix $\mathbf{C}$ contains information about the gradient's direction on
    average through its eigendecomposition. In certain cases, one observes a
    jump in the eigenvectors: a ``spectral gap." This is indicative of an
    active subspace, defined as the principal eigenspace corresponding to the
    eigenvalues prior to the jump.
	
   
   
   
    \resp{2-1}{Related to the ASM is the concept of Sufficient Dimension Reduction (SDR),
    which aims to find a subspace capturing the stochastic dependence between the independent variables
    and the response. 
    However, SDR is applied in the context of
    observational data analysis, whereas methods bearing the
    name active subspaces tend to have a UQ application in mind, where the experimenter chooses the
    inputs. Hence SDR has its roots in the multivariate statistics community.} 
	We briefly
    outline some of the most popular SDR methods and refer the reader to the
    work of \citet{yanyuan2013} for a more complete review.
	
	Sliced inverse regression \citep{li1991}, among the earliest methods,
	involves breaking the response into bins, taking the mean of each feature
	within each bin, and then performing principal component analysis on the
	bins to give important directions. \resp{}{\cite{Glaws2020} relate this technique to the ridge functions
	with which the ASM is concerned, and \cite{Li2016} discuss algorithms
	suitable to both situations when the function of interest is exactly low dimensional 
	as well as when it is only approximately so.}
	
    Arguing that eigenvectors corresponding to large eigenvalues represent
    important directions, \citet{li1992} presents a method for subspace
    dimension reduction achieved by estimating the expected Hessian $\nabla^2
    f$. Computation involves solving the generalized eigenproblem on matrices
    with size equal to the dimension of the input space. 
	
    Also typical have been methods based on estimating these or related
    quantities by using kernel smoothers. An important early article in this
    regard is that of \citet{samarov1993}, who considered estimation of
    several functions defined in terms of integrals, including the expected
    outer product of the gradient, defined exactly as $\mathbf{C}$, as well as the
    expected Hessian. In the same article, sample estimators are constructed
    for such quantities by using kernel-regression-based estimators of the
    gradient at each sample location. Practical use of these techniques is
    complicated by the need to choose a kernel bandwidth; and although methods
    exist to this end, none may be considered ``best," and the methods may
    lead to different results \citep{Ghosh2018}. \citet{Fukumizu2014}
    alleviate some of those shortcomings but still rely on a MC estimator of
    $\mathbf{C}$. Among kernel methods, GPs offer additional benefits by taking a more
    probabilistic point of view.
	
	
	\subsection{Gaussian Processes and Linear Embeddings}
	
    Gaussian Processes are a common Bayesian nonparametric technique; see, for
    example, \citet{Rasmussen2006} for an introduction. Given any set of pairs
    $(\mathbf{x}_1, y_1), \ldots, (\mathbf{x}_n, y_n)$ with $\mathbf{x}_i$ belonging to a set $\mathcal{X}$
    representing inputs and $y_i \in \mathbb{R}$ representing outputs of some
    real-valued function defined on $\mathcal{X}\subseteq \mathbb{R}^m$, the GP models
    the output vector $\mathbf{y}_n = \left(y_i \right)_{1 \leq i \leq n}$ as coming
    from a Gaussian distribution with a mean vector and covariance matrix
    dependent on $\mathbf{X}_n = (\mathbf{x}_i)_{1 \leq i \leq n}$. Typically, the mean is set
    to be a constant or a linear function of $\mathbf{X}_n$, and the covariance matrix
    is determined via a \textit{kernel function}, a positive definite function
    defined on $\mathcal{X}\times \mathcal{X}$, such that the covariance satisfies
    \resp{1-4}{$\Cov{\omega}{y_i, y_j} = k(\mathbf{x}_i, \mathbf{x}_j)$} \resp{}{(here $\omega$ is an element of the sample space $\Omega$)}. 
    GP regression (also known as kriging) may be interpreted as a Bayesian linear regression
    on an infinite-dimensional space determined by the chosen kernel function 
    \citep{Scholkopf2001}.
	
    For the purposes of this article, suppose we are given $n$ observations of
    a \resp{}{(possibly noisy)} black-box function 
\resp{}{mapping the input space $\mathcal{X}$ and the stochastic outcome space $\Omega$ to the reals} $f: \mathcal{X} \times \Omega \rightarrow \mathbb{R}$, possibly
    depending on an active subspace of dimension $r$. Based on this data \resp{}{$\mathcal{A}_n
    = \left(\mathbf{x}_i, y_i = f(\mathbf{x}_i, \omega_i)  \right)_{1 \leq i \leq n}$}, and considering a
    prior (zero-mean) Gaussian process $Y \sim GP(0, k(\cdot,\cdot))$,
    classical multivariate normal properties result in the GP predictive
    equations; that is, $Y|\mathcal{A}_n \sim \mathcal{N}(m_n(\cdot),
    k_n(\cdot,\cdot))$ with
	\begin{align}
	m_n(\mathbf{x}) &= \veck(\mathbf{x}) \mathbf{K}_n^{-1} \mathbf{y}_n\\
	k_n(\mathbf{x}, \mathbf{x}') &= k(\mathbf{x},\mathbf{x}') - \veck(\mathbf{x}) \mathbf{K}_n^{-1} \veck(\mathbf{x}')^\top,
	\label{eq:kn}
	\end{align}
	where $\veck(\mathbf{x}) = \left(k(\mathbf{x}, \mathbf{x}_i) \right)_{1 \leq i \leq n}$ and $\mathbf{K}_n
	= \left( k(\mathbf{x}_i, \mathbf{x}_j) \right)_{1 \leq i,j \leq n}.$
	
	
    Unfortunately, without modification, GPs are not generally suited for
    problems of more than a few tens of variables. A common remedy is to
    conduct some manner of dimension reduction and perform the kriging on this
    smaller space. Popular among these methods are linear embeddings, which
    involve projecting the input data onto a lower-dimensional subspace. This
    may be done randomly as by \citet{Wang2016} or may be treated as kernel
    hyperparameters. A rank-deficient matrix in GP kernels to reduce
    dimensionality is presented, e.g., by \citet{Vivarelli:1999}. That the
    matrix is rank-deficient means that its range has smaller dimension than
    the underlying space, resulting in dimension reduction.
    \citet{Tripathy2016} explicitly model the GP as existing in the
    low-dimensional space defined by $\mathbf{A}\mathbf{x}$ for some $\mathbf{A}$
    (though this is easily seen to be equivalent to putting $\mathbf{A}$ in
    the kernel) and propose a two-step approach by alternating between
    optimizing with respect to the matrix $\mathbf{A}$ and the GP kernel
    hyperparameters. Since $\mathbf{A}$ is constrained to be orthogonal for
    identifiability purposes, this involves a nontrivial optimization over
    Grassmann manifolds (see also \citep{Marcy2017}). \citet{Hokanson2018}
    describe a similar approach using a polynomial model on the active
    subspace for which the two-step algorithm simplifies. Also taking the
    hyperparameter route are \citet{Garnett2014}, who define an acquisition
    function for sequential design that minimizes uncertainty on the subspace.
    This work bears similarity to ours but differs substantially in that it
    fits a GP on a low-dimensional subspace to be estimated, whereas we fit a
    GP on the entire space and use properties of GPs to estimate the important
    subspace.
	
	\resp{2-2/3/4/5/6/7/8}{
	\subsection{Linking GPs and ASMs}
	
	GPs via linear embeddings (GPLEs) and ASMs have much in common. Although
	links between GPs and ASMs have been observed previously, here we
	establish the prior properties of the classical ASM estimator $\hat{\mathbf{C}}$ as a method of
	moments (MoM) estimator for the linear embedding matrix in GPLEs in
	Theorem \ref{th:mom} and its corollary.

    \begin{Theorem}
		Assume we are given a mean-zero Gaussian process with the stationary
		kernel
		$
		k(\mathbf{x}_i,\mathbf{x}_j) = \sigma^2 \exp\{-0.5 (\mathbf{x}_i-\mathbf{x}_j)^\top\mathbf{A}^\top\mathbf{A}(\mathbf{x}_i-\mathbf{x}_j)\}
		$,
		that is, an anisotropic Gaussian kernel with variance parameter
		$\sigma^2$ and Mahalanobis distance given by the matrix
		$\mathbf{A}^\top\mathbf{A}$, with $\mathbf{A} \in
		\mathbb{R}^{r\times m}$.
		Let $f \sim \textrm{GP}(0, k(\cdot, \cdot))$. 
		Then,
		$
		\Esp{\omega}{(\nabla f(\mathbf{x}_i, \omega)) (\nabla f(\mathbf{x}_i, \omega))^\top} = \sigma^2 \mathbf{A}^\top\mathbf{A}
		$
		for any $\mathbf{x}_i$, where $\omega$ denotes the source of stochasticity in $f$.
		\label{th:mom}
    \end{Theorem}
	\begin{proof}
	
    Since the active subspace is defined on the gradients of the function of
    interest, it is useful to derive the GP implied by the above model on the
    gradient.
		
    That the stochastic process implied by a GP on a function's partial
    derivatives is also a GP is well known.
    In particular, the prior covariance implied on partial derivatives of the GP is given by second-order derivatives
    of the kernel function \citep[Chapter 9.4]{Rasmussen2006}:
	\begin{equation*}
		\Cov{\omega}{\frac{\partial f(\mathbf{x},\omega)}{\partial x_{i,d}}, 
		\frac{\partial f(\mathbf{x},\omega)}{\partial x_{j,e}}} =
		\frac{\partial^2 k(\mathbf{x}_{i}, \mathbf{x}_{j})}{\partial x_{i,d}\partial x_{j,e}}.
	\end{equation*}

	where $e,d\in\{1,\ldots,m\}$ give input variables.
	It is therefore useful to derive the Hessian of our kernel function. Given
	that the derivative is
	\begin{equation*}
		\frac{\partial k(\mathbf{x}_1, \mathbf{x}_2)}{\partial \mathbf{x}_1} = 
		- \sigma^2 \mathbf{A}^\top\mathbf{A}(\mathbf{x}_1 - \mathbf{x}_2)  
		\exp\left\{-0.5 (\mathbf{x}_1 - \mathbf{x}_2)^\top\mathbf{A}^\top\mathbf{A}(\mathbf{x}_1 - \mathbf{x}_2)\right\},
	\end{equation*}		
	the Hessian can be expressed as
	\begin{equation*}
			\frac{\partial^2 k(\mathbf{x}_1, \mathbf{x}_2)}{\partial \mathbf{x}_1\partial 
\mathbf{x}_2^\top} = 
		\sigma^2 \left(\mathbf{I} - \mathbf{A}^\top\mathbf{A}(\mathbf{x}_1 - 
\mathbf{x}_2)(\mathbf{x}_1 - \mathbf{x}_2)^\top\right)
		\mathbf{A}^\top\mathbf{A}
		\exp\left\{-0.5 (\mathbf{x}_1 - \mathbf{x}_2)^\top\mathbf{A}^\top\mathbf{A}(\mathbf{x}_1 - 
\mathbf{x}_2)\right\}.
	\end{equation*}
		
	Now turn to the quantity of interest, 
	\begin{equation*}
		\Esp{\omega}{\nabla f(\mathbf{x},\omega)\nabla f(\mathbf{x},,\omega)^\top} =
		\Var{\omega}{\nabla f(\mathbf{x},\omega)} + \Esp{\omega}{\nabla f(\mathbf{x},\omega)}\Esp{\omega}{\nabla f(\mathbf{x},\omega)}^\top.
	\end{equation*}
	
	Using our zero-mean assumption for our GP, our gradient's mean is then
	determined to be zero as well
	\resp{2-6}{\footnote{$\Esp{\omega}{\frac{\partial f(\mathbf{x},\omega)}{\partial x_{d}}}=
			\frac{\partial \Esp{\omega}{f(\mathbf{x},\omega)}}{\partial x_{d}} =
			\frac{\partial 0}{\partial x_{d}} =
			0$}}.
	 This deletes the second term, leaving
	\begin{equation*}
		\Esp{\omega}{\nabla f(\mathbf{x},\omega)\nabla f(\mathbf{x},\omega)^\top} =
		\Var{\omega}{\nabla f(\mathbf{x},\omega)}. 
	\end{equation*}
	
	which we can evaluate using the Hessian derived above:
	
	\begin{align*}
		\Var{\omega}{\nabla f(\mathbf{x},\omega)} &= \Cov{\omega}{\nabla f(\mathbf{x},\omega), \nabla f(\mathbf{x},\omega)}  \\
		&= \sigma^2 (\mathbf{I} - \mathbf{A}^\top\mathbf{A}(\mathbf{x} - \mathbf{x})(\mathbf{x} - \mathbf{x})^\top)
		\mathbf{A}^\top\mathbf{A}
		\exp\{-0.5 (\mathbf{x} - \mathbf{x})^\top\mathbf{A}^\top\mathbf{A}(\mathbf{x} - \mathbf{x})\} \\
		&= \sigma^2 \mathbf{A}^\top\mathbf{A}.
	\end{align*}

	\end{proof}

	\begin{Corollary}
	For any sample $\{\mathbf{x}_i\}$, $\Esp{\omega}{\frac{1}{n} \sum_{i =1 }^{n} \nabla f(\mathbf{x}_i,\omega) \nabla f(\mathbf{x}_i,\omega)^\top} =  \sigma^2 \mathbf{A}^\top\mathbf{A}$.
	\end{Corollary}
	
	Note that no assumption was made on the measure with respect to which the inputs $\mathbf{x}$ were sampled. This is a consequence of our stationary assumption: we are effectively assuming that the directions of global importance are important at each input point. To do otherwise would require a nonstationary kernel. This section has covered estimation of GP hyperparameters via the active subspace estimator $\hat{\mathbf{C}}$. There have also been efforts to, somewhat conversely, estimate $\mathbf{C}$ from a standard Gaussian process without linear embeddings. We will review this in the next section and further develop it in the remainder of the article.

}
	\subsection{Surrogate-Assisted Active Subspace Estimation}
	
    Estimating the matrix $\mathbf{C}$ is not straightforward if gradients for the
    function $f$ are not available, as is often the case in practice. Along
    with finite-difference gradient estimation, several gradient-free methods
    have been proposed to deal with this issue given $n$ observations of
    design points stored in a matrix $\mathbf{X}_n$ and the function $f$ evaluated at
    those points stored in a vector $\mathbf{y}_n$.
	
    Perhaps the simplest of these ideas is proposed by
    \citet{Constantine2015}: simply perform linear regression of $\mathbf{y}_n$ upon
    $\mathbf{X}_n$ and look at the subspace generated by the estimated regression
    coefficient vector. This method, applicable only when a one-dimensional
    subspace is desired, was found to be adequate experimentally on several
    problems. Also proposed by \citet{Constantine2015} is the use of local
    linear models, together with methods for aggregating their coefficient
    vectors. \resp{2-9}{ Further, \cite{constantine2015sketch} show how to estimate the active subspace 
   	when the full gradient is unavailable but certain directional derivatives are known.}
	
    Another approach is to retain the MC estimator but to use surrogate
    modeling to estimate the gradient at that point. \citet{palar2017},
    \citet{Othmer2016} and \citet{Namura2017} use finite differencing on a GP
    or polynomial chaos expansion in order to estimate gradients at random
    points throughout the input space, using the MC estimate of $\mathbf{C}$.
	
    We now show that the GP assumptions make such a process unnecessary, since
    a closed-form estimate is available.
	

	\section{Methodology}
	\label{sec:methodology}
	
	We now derive a closed-form estimator and a novel methodology for its use
	in sequential procedures. Our analytic estimator eases the computational
	burden and gives room for other procedures to be built on it.
	
	\subsection{Closed-Form Estimator of \texorpdfstring{$\mathbf{C}$}{C} from GPs}
	\label{sec:CGP}
    Recall our fundamental inference task of estimating $\mathbf{C} =  \Esp{\mathbf{x}}{(\nabla f)
    (\nabla f)^\top}$ given a set of observations $(\mathbf{X}_n, \mathbf{y}_n)$. Although
    $\mu$ will be the Lebesgue measure in our experiments (Section
    \ref{sec:experiments}), we will neither make nor need the assumption that
    $\mathbf{X}_n \sim \mu$. Given our Gaussian assumptions on the stochastic process,
    the expectation is available analytically in terms of standard functions
    for popular kernels.

    Let us now express the matrix $\C^{(n)}$ for a GP, starting \resp{2-10}{with 
    design points 
    $\mathcal{A}_n = \left(\mathbf{x}_i, y_i = f(\mathbf{x})  \right)_{1 \leq i \leq n}$
	}
     before providing sequential expressions.
    Recall that $\mathbf{K}_n$ is the kernel matrix given $n$ observations: given a
    positive definite kernel function $k$, the components of $\mathbf{K}_n$ are given
    by $(\mathbf{K}_n)_{i,j} = k(\mathbf{x}_i, \mathbf{x}_j) + \tau^2 \delta_{i=j}$, where, in order
    to also account for noisy observations, $\delta$ is the Kronecker delta
    and $\tau^2$ is the error variance.
	
    Differentiation being a linear operator, it is well known that $\nabla Y =
    \partial Y / \partial x_1, \dots, \partial Y / \partial x_m$ is also a
    Gaussian process. Assuming that k is twice differentiable, the joint
    distribution of $\left(Y(\mathbf{X}_n), \partial Y(\mathbf{x}) / \partial \mathbf{x}_1,
    \dots, \partial Y(\mathbf{x}) / \partial \mathbf{x}_m \right)$ is 
	\[
	\begin{pmatrix}
	\mathbf{y}_n\\
	\partial Y(\mathbf{x}) / \partial \mathbf{x}_1\\
	\vdots\\
	\partial Y(\mathbf{x}) / \partial \mathbf{x}_m
	\end{pmatrix}
	\sim \mathcal{N}
	\left(
	\begin{pmatrix}
	\mathbf{0_n}\\
	0\\
	\vdots\\
	0
	\end{pmatrix},
	\begin{pmatrix}
	\mathbf{K}_n & \partial \veck(\mathbf{x})^\top / \partial x_1 & \hdots & \partial \veck(\mathbf{x})^\top / \partial x_m\\
	\partial \veck(\mathbf{x}) / \partial x_1 & \partial^2k(\mathbf{x}, \mathbf{x})/\partial x_1^2 & \hdots &  \partial^2k(\mathbf{x}, \mathbf{x})/\partial x_1 \partial x_m\\
	\vdots & \vdots & \ddots & \vdots\\
	\partial \veck(\mathbf{x}) \partial / x_m & \partial^2k(\mathbf{x}, \mathbf{x})/\partial x_m \partial x_1 & \hdots & \partial^2 k(\mathbf{x}, \mathbf{x})/\partial x_m^2\\
	\end{pmatrix}
	\right),
	\]
	or, in shorthand,
	\[
	\begin{pmatrix}
	\mathbf{y}_n\\
	\nabla Y (\mathbf{x})
	\end{pmatrix}
	\sim \mathcal{N}
	\left(
	\begin{pmatrix}
	\mathbf{0_n}\\
	\mathbf{0_m}
	\end{pmatrix},
	\begin{pmatrix}
	\mathbf{K}_n & \boldsymbol{\kappa}(\mathbf{x})^\top\\
	\boldsymbol{\kappa}(\mathbf{x}) & \mathbf{K}_m(\mathbf{x},\mathbf{x}) 
	\end{pmatrix}
	\right),
	\]
	wherein $\boldsymbol{\kappa}_i(\mathbf{x}) = \frac{\partial \veck(\mathbf{x})}{\partial x_i}\in\mathbb{R}^n$, for $i=1, \ldots, m$.


	As a result, $\nabla Y(\mathbf{x})|\mathcal{A}_n \sim \mathcal{N}(\boldsymbol{\mu}_n(\mathbf{x}), \tilde{\mathbf{K}}_n(\mathbf{x}, \mathbf{x}))$ with 
	\begin{align*}
	\boldsymbol{\mu}_n(\mathbf{x}) &= \boldsymbol{\kappa}(\mathbf{x}) \mathbf{K}_n^{-1} \mathbf{y}_n\\
	\tilde{\mathbf{K}}_n(\mathbf{x}, \mathbf{x}) & = \mathbf{K}_m(\mathbf{x}, \mathbf{x}) - \boldsymbol{\kappa}(\mathbf{x}) \mathbf{K}_n^{-1} 
	\boldsymbol{\kappa}(\mathbf{x})^\top,
	\end{align*} 
	which leads to a closed-form expression for $\mathbf{C}$, as given in Theorem \ref{th:C_GP}.
	\begin{Theorem}
        Let $k$ be a twice differentiable kernel, $\mathbf{W}_{i,j} :=
		\int_\mathcal{X} \boldsymbol{\kappa}_i(X)^\top \boldsymbol{\kappa}_j(X) d\mu$,
		and $E_{i,j} :=
		\int_\mathcal{X} \frac{\partial^2 k (X, X)}{\partial x_i \partial x_j}
		d\mu$. Then,  
		$C^{(n)}_{i,j} = E_{i,j} - tr\left( \mathbf{K}_n^{-1} \mathbf{W}_{i,j} \right) +
		\mathbf{y}_n^\top \mathbf{K}_n^{-1} \mathbf{W}_{i,j} \mathbf{K}_n^{-1} \mathbf{y}_n$.
		\label{th:C_GP}
	\end{Theorem}

	\begin{proof}
		Based on \citet[Equation (321)]{Petersen2008}, $\C^{(n)} = \Esp{\mathbf{x}, \omega}{\nabla
		Y(X)(\nabla Y(X))^\top |\mathcal{A}_n } = \mathbf{M} + \mathbf{m} \mathbf{m}^\top$ with $X$ a
		random vector (e.g., same distribution as used for sampling the
		gradient),  $\mathbf{m} = \Esp{\mathbf{x},\omega}{\nabla Y(X)|\mathcal{A}_n }$, $\mathbf{M} = \Var{\mathbf{x},\omega}{\nabla
		Y(X)|\mathcal{A}_n }$. We use integration and expectations
		interchangeably.
		
		The first term is 
		$\mathbf{m} = \Esp{\mathbf{x}}{\boldsymbol{\mu}_n(X)} = \Esp{\mathbf{x}}{\boldsymbol{\kappa}(X)} \mathbf{K}_n^{-1} \mathbf{y}_n 
		= \left[\int_{\mathcal{X}} \boldsymbol{\kappa}(\mathbf{x}) d\mu \right] \mathbf{K}_n^{-1} \mathbf{y}_n$.
		For the second one, $\mathbf{M}$, we use the law of total covariance:
		\begin{align*}
		M_{i,j} =& \Cov{\mathbf{x},\omega}{ \frac{\partial Y(X)}{\partial x_i}, 
		\frac{\partial Y(X)}{\partial x_j} |\mathcal{A}_n}\\
		=& \Esp{\mathbf{x}}{\Cov{\omega}{ \frac{\partial Y(X)}{\partial x_i}, 
		\frac{\partial Y(X)}{\partial x_j} |\mathcal{A}_n, X } } + \Cov{\mathbf{x}}{\Esp{\omega}{ 
		\frac{\partial Y(X)}{\partial x_i} | \mathcal{A}_n, X}, 
		\Esp{\omega}{\frac{\partial Y(X)}{\partial x_j} | \mathcal{A}_n, X} } \\ 
		=& \Esp{\mathbf{x}}{\tilde{K}_{i,j}^{(n)}(X,X)}+ \Cov{\mathbf{x}}{ \mu_{i}^{(n)}(X), \mu_{j}^{(n)}(X) }\\
		=&  \Esp{\mathbf{x}}{\frac{\partial^2 k (X, X)}{\partial x_i \partial x_j} - 
		\boldsymbol{\kappa}_i(X) \mathbf{K}_n^{-1} \boldsymbol{\kappa}_j(X)^\top} + 
		\Cov{\mathbf{x}}{\boldsymbol{\kappa}_i(X) \mathbf{K}_n^{-1} \mathbf{y}_n,
		\boldsymbol{\kappa}_j(X) \mathbf{K}_n^{-1} \mathbf{y}_n } \\
		=& E_{i,j} - \Esp{\mathbf{x}}{\boldsymbol{\kappa}_i(X) \mathbf{K}_n^{-1} \boldsymbol{\kappa}_j(X)^\top} 
		+ \mathbf{y}_n^\top \mathbf{K}_n^{-1} \Cov{\mathbf{x}}{\boldsymbol{\kappa}_i(X),
		 \boldsymbol{\kappa}_j(X)}\mathbf{K}_n^{-1} \mathbf{y}_n\\
		=& E_{i,j} - tr \left(\mathbf{K}_n^{-1} \left[\Cov{\mathbf{x}}{ \boldsymbol{\kappa}_i(X), 
		\boldsymbol{\kappa}_j(X) } + \Esp{\mathbf{x}}{\boldsymbol{\kappa}_i(X)}^\top 
		\Esp{\mathbf{x}}{\boldsymbol{\kappa}_j(X)} \right] \right) +\\
		& \mathbf{y}_n^\top \mathbf{K}_n^{-1} \Cov{\mathbf{x}}{\boldsymbol{\kappa}_i(X), 
		\boldsymbol{\kappa}_j(X)} \mathbf{K}_n^{-1} \mathbf{y}_n,
		\end{align*}
		where 
		$\boldsymbol{\kappa}_i(\mathbf{x}) = \Cov{\mathbf{x}}{ Y(\mathbf{X}_n), \frac{\partial Y (\mathbf{x})}{\partial x_i}}$.
		
		The result follows by writing $C^{(n)}_{i,j} = M_{i,j} + m_i m_j$ and 
		denoting $\mathbf{W}_{i,j} = \Cov{\mathbf{x}}{ \boldsymbol{\kappa}_i(X), \boldsymbol{\kappa}_j(X)} + 
		\Esp{\mathbf{x}}{\boldsymbol{\kappa}_i(X)}^\top \Esp{\mathbf{x}}{\boldsymbol{\kappa}_j(X)} = \int_\mathcal{X} \boldsymbol{\kappa}_i(X) \boldsymbol{\kappa}_j(X)^\top d\mu$.
	\end{proof}
	
    Hence, once a form for the kernel is chosen \resp{2-12}{and the data have been observed,}
     the active subspace of the GP
    depends only on the kernel hyperparameters (e.g., the $m$
    length scales), in contrast with the matrix parameterizations of
    \cite{Garnett2014} and \cite{Tripathy2016} that involve $m \times r$
    hyperparameters. For classically used kernels such as tensor product
    versions of the Gaussian and Mat\'ern (with smoothness parameter $\nu \in
    \{3/2, 5/2\}$) kernels, $E_{i,j}$ and $\mathbf{W}_{i,j}$ are available even in closed
    form (see Appendix \ref{ap:kernelexps} and supplementary materials).
    $C^{(n)}_{i,j}$ can be decomposed into two parts: the first two terms
    correspond to an integrated mean square prediction error (especially
    apparent when $i = j$), while the second is related to the product of
    predictive means. \resp{2-12}{The matrix $\C^{(n)}$ is thus estimated via 
    the (nonstationary, in contrast to Theorem \ref{th:mom}) posterior Gaussian process
    conditioned on the observed input-output pairs $\mathcal{A}_n$. }
    	
   
   
    
    
    The diagonal terms of the matrix can also be connected
    to derivative-based global sensitivity measures, for which the expressions
    have also been derived by \citet{DeLozzo2016}.

    Given our focus here on sequential design, also of interest is an update
    formula for calculating $\C^{(n+1)}$ (i.e., our estimate for $\mathbf{C}$ given
    $\mathcal{A}_{n+1}$) accounting for the fact that we have already calculated
    $\C^{(n)}$ using $\mathbf{K}_n$ and $\mathbf{W}_{i,j}^{(n)}$. Specifically, $\mathbf{K}_{n+1}$ and $\mathbf{W}_{i,j}^{(n+1)}$
    can be decomposed as
    \[
	\mathbf{K}_{n+1}=\begin{bmatrix}
	\mathbf{K}_{n} & \veckn(\tilde{\x}) \\
	\veckn(\tilde{\x})^\top & k(\tilde{\x},\tilde{\x}) + \tau^2
	\end{bmatrix}, \quad
	\mathbf{W}_{i,j}^{(n+1)}= \begin{bmatrix}
	\mathbf{W}_{i,j}^{(n)} & \mathbf{w}_a(\tilde{\x})\\
	\mathbf{w}_b(\tilde{\x})^\top & w(\tilde{\x}, \tilde{\x})
	\end{bmatrix}.
	\] 

    The resulting update formula is given in Theorem \ref{th:Cn1}.
	\begin{Theorem}
	Denote $\mathbf{g}(\tilde{\x})= -\sigma_n^2(\tilde{\x})^{-1}\mathbf{K}_n^{-1} \veckn(\tilde{\x})$ and
	$\sigma^2_n(\tilde{\x}) = k_n(\tilde{\x}, \tilde{\x})$. Given a new design point $\tilde{\x}$
	but not the function value at this location (i.e., $y_{n+1} \sim
	\mathcal{N}(m_n(\tilde{\x}), k_n(\tilde{\x}, \tilde{\x}))$), the random variable
	$C_{i,j}^{(n+1)} - C_{i,j}^{(n)}$  can be written as
	\begin{equation*}
	 C_{i,j}^{(n+1)} - C_{i,j}^{(n)} = \alpha_{i,j}(\tilde{\x}) + 
	 Z \beta_{i,j}(\tilde{\x}) + Z^2 \gamma_{i,j}(\tilde{\x}),
	\end{equation*}
	with $Z \sim \mathcal{N}(0,1)$ and
		\begin{align*}
	\alpha_{i,j}(\tilde{\x}) =& - (\mathbf{w}_a(\tilde{\x}) + \mathbf{w}_b(\tilde{\x}))^\top \mathbf{g}(\tilde{\x}) -
	 \frac{w(\tilde{\x}, \tilde{\x}) + 
	 \veckn(\tilde{\x})^\top \mathbf{K}_n^{-1} \mathbf{W}_{i,j}^{(n)} \mathbf{K}_n^{-1} \veckn(\tilde{\x})}{\sigma^2_n(\tilde{\x})}\\ 
	\beta_{i,j}(\tilde{\x}) =&  \frac{\mathbf{y}_n^\top \mathbf{K}_n^{-1} 
	\mathbf{W}_{i,j}^{(n)} \mathbf{K}_n^{-1} \veckn(\tilde{\x}) + \veckn(\tilde{\x})^\top \mathbf{K}_n^{-1} \mathbf{W}_{i,j}^{(n)} \mathbf{K}_n^{-1} 
	\mathbf{y}_n - (\mathbf{w}_a(\tilde{\x}) + \mathbf{w}_b(\tilde{\x}))^\top \mathbf{K}_{n}^{-1} \mathbf{y}_n}{\sigma_n(\tilde{\x})} \\ 
	\gamma_{i,j}(\tilde{\x}) =&   \frac{w(\tilde{\x}, \tilde{\x}) 
	+ \veckn(\tilde{\x})^\top \mathbf{K}_n^{-1} \mathbf{W}_{i,j}^{(n)} \mathbf{K}_n^{-1} \veckn(\tilde{\x})- (\mathbf{w}_a(\tilde{\x}) + 
	\mathbf{w}_b(\tilde{\x}))^\top \mathbf{K}_n^{-1} \veckn(\tilde{\x})}{\sigma_n^2(\tilde{\x})}.
	\end{align*}
	Once the response $\mathbf{y}_{n+1}$ has been observed, the expression evaluates to
	\begin{align*}
	&C_{i,j}^{(n+1)} - C_{i,j}^{(n)}= - (\mathbf{w}_a(\tilde{\x}) + \mathbf{w}_b(\tilde{\x}))^\top \mathbf{g}(\tilde{\x})\\
	&- \left(\mathbf{y}_n^\top \mathbf{g}(\tilde{\x}) + \frac{y_{n+1}}{\sigma_n^2(\tilde{\x})}\right) \left[ \mathbf{y}_n^\top \mathbf{K}_n^{-1} \mathbf{W}_{i,j}^{(n)} \mathbf{K}_n^{-1} \veckn(\tilde{\x}) + \veckn(\tilde{\x})^\top \mathbf{K}_n^{-1} \mathbf{W}_{i,j}^{(n)} \mathbf{K}_n^{-1} \mathbf{y}_n \right]\\ 
	&+ \left(\mathbf{y}_n^\top \mathbf{g}(\tilde{\x}) + \frac{y_{n+1}}{\sigma_n^2(\tilde{\x})}\right) (\mathbf{w}_a(\tilde{\x}) + \mathbf{w}_b(\tilde{\x}))^\top (\mathbf{K}_{n}^{-1} \mathbf{y}_n + \mathbf{g}(\tilde{\x}) y_{n+1} - \mathbf{g}(\tilde{\x}) \veckn(\tilde{\x})^\top \mathbf{K}_n^{-1} \mathbf{y}_n)\\
	&+ \left[ \left(\mathbf{y}_n^\top \mathbf{g}(\tilde{\x}) + \frac{y_{n+1}}{\sigma_n^2(\tilde{\x})}\right)^2 - 
	\frac{1}{\sigma_n^2(\tilde{\x})} \right] \left[ w(\tilde{\x}, \tilde{\x}) + \veckn(\tilde{\x})^\top 
	\mathbf{K}_n^{-1} \mathbf{W}_{i,j}^{(n)} \mathbf{K}_n^{-1} \veckn(\tilde{\x}) \right].
	\end{align*}
	\label{th:Cn1}
	\end{Theorem}
	\begin{proof}
	The derivation is detailed in Appendix \ref{sec:update}.
	\end{proof}

	Note that in contrast to the classical MC estimator $\hat{\mathbf{C}}$, by using a
	GP to predict everywhere, one can follow an ``off-policy" approach: design
	points $\mathbf{X}_n$ need not be sampled from $\mu$. Combined with the above update
	expression, this enables sequential design capabilities to learn active
	subspaces more efficiently.
	
		
	\subsection{Illustration: Comparison with ARD}
	\label{ard}
	
    At first glance, a disadvantage of kernel methods as compared with
    classical parametric models is the difficulty in interpreting model
    parameters. A popular heuristic for determining relative variable
    importance is that of automatic relevance determination (ARD)
    \citep[Chapter 5.1]{Rasmussen2006}, which asserts that input dimensions
    with large length scales are of little importance. This is due to the
    limiting behavior of many popular kernels that asymptotically ignore an
    input dimension as its length scale tends to infinity. This is taken
    advantage of, for example, for the purpose of axis-aligned dimension
    reduction for optimization \citep{Salem2018}. For two finite length scales
    $l_1 > l_2$ belonging to input variables $x_1, x_2$, however, there is no
    guarantee that $x_2$ is more important than $x_1$. As we illustrate here
    with a simple example, examining the loadings of variables in the active
    subspace eigenbasis may be more informative.
	
    Consider the function 
    \begin{equation}
     f(x_1, x_2) = a \sin(bx_1) + c x_2^2; \qquad a = 0.1, \, b = 20, \, c = -4;
     \qquad (x_1,x_2) \in [0,1]^2
    \label{eq:testfun}
     \end{equation}
    illustrated in Figure~\ref{fig:ard}. Along the $x_1$ dimension is a
    sinusoidal function with high frequency but low amplitude, and along the
    other is simply a quadratic function in $x_2$, which dominates in terms of
    determining the function value. Given 1,000 observations uniformly
    distributed from the 2D unit interval, a GP with a Gaussian kernel gives
    length-scale estimates (determined via maximum likelihood estimation) of
    0.069 for $x_1$ and 0.37 for $x_2$, which, according to the ARD principle,
    mistakenly suggests $x_1$ is more important than $x_2$.
	
    Turning now to the active subspace estimated by the GP for the same
    sample, we find the eigendecomposition of our estimated $\mathbf{C}$ to be
	\begin{equation*}
		\hat{\mathbf{C}} = \begin{bmatrix}
		-0.00 & -0.99 \\
		0.99 & -0.00
		\end{bmatrix}
		\quad
		\begin{bmatrix}
		13.63 & 0 \\
		0 & 1.30
		\end{bmatrix}
		\quad
		\begin{bmatrix}
		-0.00 & 0.99 \\
		-0.99 & -0.00
		\end{bmatrix}.
	\end{equation*}
    The first eigenvector represents $x_2$ while the second represents $x_1$,
    and the first eigenvalue is an order of magnitude greater than the first,
    suggesting, correctly, that the variable $x_2$ is of greater importance
    \footnote{\resp{2-13}{A similar example was considered in \cite{constantine2017stationary}, where the magnitudes of $a$ and $c$ were switched, and the opposite conclusion of variable importance was reached.}}.
    The effect of projecting along the eigenvectors is further illustrated in
    Figure~\ref{fig:ard}. 
    \resp{1-6}{The ASM is not the only global sensitivity
    metric that chooses the ``right'' coordinate: indeed, the Sobol indices \citep{sobol2001} for our two variables are 
    approximately $0.003$ and $0.997$, respectively. 
    The power in the ASM lies in its ability to detect non-axis-aligned directions of importance.
    }
	
	\begin{figure}[ht]
		\centering
		\includegraphics[width=\textwidth]{images/counter_proj.pdf}
		\caption{The 2-D function \resp{1-7}{given by Equation} \eqref{eq:testfun} 
		exhibits insignificant but high-frequency behavior in one dimension 
		and significant but long-range behavior in the other. Filled black
		points are design points while red circles represent points evaluated
		on a regular grid for visualization. The ASM
		reveals a linear subspace along which most of the important changes
		occur (center), while the ARD principle would select the
		high-frequency dimension (right).}
        \label{fig:ard}
	\end{figure}
		
	\subsection{Quantifying the Uncertainty in C}
	
    Of perhaps equal importance as giving an estimate for $\mathbf{C}$ is quantifying
    uncertainty in that estimate. Given an estimate, \citet{Constantine2015}
    shows how the bootstrap method \citep{Efron1981} can be used to create
    intervals by applying parametric bootstrapping to simple methods used to
    estimate active subspaces (e.g., OLS, LL).
    However, application of the bootstrap to more sophisticated active
    subspace estimation methods such as the one presented in this article can
    lead to computational challenges. Instead, uncertainty of GP parameters
    can be propagated to uncertainty of $\mathbf{C}$ via a simple MC procedure. Myriad
    methods 
    for quantifying GP hyperparameter uncertainty exist: a full
    posterior can be developed by using MCMC procedures, approximate
    posteriors can be determined (e.g., via variational methods), or the
    asymptotic normality of maximum likelihood estimators can be exploited. In
    light of its computational ease, the last method is explored here. This
    method involves computing the Hessian of the log likelihood at its maximum
    value, giving a covariance matrix for the parameters, and treating the
    maximizing parameter setting as the mean. This uncertainty on model
    parameters may be propagated to the active-subspace-defining singular
    values via MC: for each drawn parameter setting, estimate the singular
    values, then keep, for instance, the middle 95 percent of the draws. As an
    illustration, in Figure~\ref{fig:uq} we form intervals on the eigenvalues
    of the problem defined in \eqref{eq:testfun}. The downside of this
    approach is that, while it is speedy, it does not account for uncertainty
    as well as would a fully Bayesian approach leveraging MCMC.
    
    \begin{figure}[htpb]
    	\includegraphics[scale=0.45]{images/eig_ci_20.pdf}
    	\includegraphics[scale=0.45]{images/eig_ci_40.pdf}
	   	\caption{Eigenvalues from the problem in \eqref{eq:testfun}
	   	are given interval estimates, expressing decreasing model uncertainty
    	with samples of size 20 (left) and 40 (right). The Hessian expansion
    	of the log-likelihood at its maximum gives a correlation of 0.16
    	between the two length-scale parameters. 95 and 99\% confidence
    	intervals are given in black and red, respectively, while the
    	horizontal lines give the locations of the ``true" eigenvalues
    	estimated via classical MC techniques with 10,000 samples.}
	  	\label{fig:uq}
    \end{figure}
	
	\subsection{Acquisition Functions for Active Subspaces}
	
    Equipped with a closed-form GP estimate of $\mathbf{C}$ for a fixed design, we
    turn to the problem of choosing the next design to be evaluated given our
    observed responses so as to best learn about the active subspace. As do
    \citet{Labopin-Richard2016}, we take the approach of selecting the design
    that maximizes variance of our unknown quantity $\mathbf{C}$ at the next
    iteration. However, how to define this variance in our case is not clear.
    One could define subspace variance through ideas such as Fr\'echet
    variance \citep{frechet1948} and seek to minimize this, but the dimension
    of the active subspace is not known ahead of time, which would have to be
    determined by human observer or via a heuristic at each step of the
    sequential design process. Also possible is to consider measuring the
    variance of the spectral gap, which is problematic for the same reasons.
    We instead focus on the variance of the matrix $\C^{(n+1)}$, which we measure
    in three ways.
	
	Since examination of $\mathbf{C}$'s spectrum is central to making decisions
	regarding the dimension of the active subspace, pinning down $\mathbf{C}$'s
	eigenvalues is of great interest. As such, we define the \texttt{Trace}
	acquisition function using the variance of the mean eigenvalue, which may be
	optimized by maximizing the variance of the trace of $\C^{(n+1)}$:
	\begin{equation*}
		\texttt{Trace} = \Var{\resp{}{\omega}}{\textrm{tr}(\C^{(n+1)})}.
	\end{equation*}
	
	Of practical concern is that the \texttt{Trace} function may be calculated
	without computing the off-diagonal elements of the $\mathbf{C}$ matrix, which
	seems to some degree antithetical to the principle of ASMs seeking
	non-axis-aligned directions of importance.
	

	
	We also consider two acquisition functions attempting to directly measure
	$\mathbf{C}$'s variance, which we term \texttt{Var1},
	\begin{equation*}
		\texttt{Var1} = \left\|\Esp{\resp{}{\omega}}{\left(\C^{(n+1)} - \Esp{\resp{}{\omega}}{\C^{(n+1)}}\right)\odot\left(\C^{(n+1)} - 
		\Esp{\resp{}{\omega}}{\C^{(n+1)}}\right)}\right\|_F^2,
	\end{equation*}
	 and \texttt{Var2},
	\begin{equation*}
		\texttt{Var2} = \left\|\Esp{\resp{}{\omega}}{\left(\C^{(n+1)} - \Esp{\resp{}{\omega}}{\C^{(n+1)}}\right)\left(\C^{(n+1)} - 
		\Esp{\resp{}{\omega}}{\C^{(n+1)}}\right)}\right\|_F^2.
	\end{equation*}
    Here, $\odot$ represents elementwise (Hadamard product) multiplication and
    $\|\mathbf{A}\|_{F}$ denotes the Frobenius norm of $\mathbf{A}$ (the root
    of the sum of its squared elements).

    Of practical interest is that, again, closed-form expressions can be
    obtained for these three acquisition functions.

	\begin{Theorem}
	Given the coefficients $\beta_{i,j}(\tilde{\x})$ and $\gamma_{i,j}(\tilde{\x})$ from
	Theorem~\ref{th:Cn1} stored in matrices $\mathbf{B}(\tilde{\x})$ and
	$\bm{\Gamma}(\tilde{\x})$, the acquisition functions and their gradients are
	available in the closed forms shown in Table~\ref{tab:acq_expr}.
	\label{th:anltAcq}
	\end{Theorem}
	\begin{proof}
		The derivation is detailed in Appendix~\ref{sec:infill}. 
	\end{proof}

	We note further that $\mathbf{B}(\tilde{\x})$ and $\bm{\Gamma}(\tilde{\x})$ have
	closed-form expressions for popular kernels, including the Gaussian,
	Mat\'ern 3/2, and Mat\'ern 5/2; see Appendix \ref{ap:kernelexps}.

	\begin{table}[h]
		\centering
		
		\caption{Analytic expressions and gradients with respect to the new
		location $\tilde{\x}$ for each of our acquisition functions. For
		compactness, we employ $\mathbf{A}}%^{(n+1)}=\C^{(n+1)} - \Esp{\resp{}{\omega}}{\C^{(n+1)}}$. Dependence of
		$\mathbf{B}(\tilde{\x}), \partial\mathbf{B}(\tilde{\x}), \bm{\Gamma}(\tilde{\x})$,
		and $\partial\bm{\Gamma}(\tilde{\x})$ on $\tilde{\x}$ has been suppressed for
		brevity. \label{tab:acq_expr}}
	


	\vspace{1em}

		\resizebox{\columnwidth}{!}{
			\begin{tabular}{l*{3}{c}}
				\textbf{Name}              & \textbf{Definition} & 
				\textbf{Expression} & \textbf{d-Derivative}\\
				\hline
				\texttt{Trace} & $\Var{\resp{}{\omega}}{\textrm{tr}(\mathbf{C}^{(n+1)})}$  & 
				$\textrm{tr}(\mathbf{B}\odot\mathbf{B}) + 
				2\textrm{tr}(\bm{\Gamma}\odot\bm{\Gamma})$ & 
				$2\textrm{tr}(\partial\mathbf{B}_d)\textrm{tr}(\mathbf{B}) + 
				4\textrm{tr}(\partial\bm{\Gamma}_d\textrm{tr}(\bm{\Gamma}))$\\
				\texttt{Var1}            & $\|\Esp{\resp{}{\omega}}{\mathbf{A}}%^{(n+1)}\odot\mathbf{A}}%^{(n+1)}}\|_F^2$ & 
				$\|\mathbf{B}\odot\mathbf{B}+ 2 \bm{\Gamma}\odot\bm{\Gamma}\|_F^2$ &  
				$2(2\mathbf{B}\odot\partial\mathbf{B}_d + 4 
				\bm{\Gamma}\odot\partial\bm{\Gamma}_d) \|\mathbf{B}\odot\mathbf{B}+ 2 
				\bm{\Gamma}\odot\bm{\Gamma}|\|_F^2$\\
				\texttt{Var2}           & $\|\Esp{\resp{}{\omega}}{\mathbf{A}}%^{(n+1)}\An}\|_F^2$ & 
				$\|\mathbf{B}\mathbf{B} + 2 
				\bm{\Gamma}\bm{\Gamma}\|_F^2$ & $2(\mathbf{B}\partial\mathbf{B}_d + 
				\partial\mathbf{B}_d\mathbf{B} + 2\bm{\Gamma}\partial\bm{\Gamma}_d + 
				2\partial\bm{\Gamma}_d\bm{\Gamma})\|\mathbf{B}\mathbf{B} + 2 
				\bm{\Gamma}\bm{\Gamma}\|_F^2$\\
			\end{tabular}
		}
	\end{table}


    The details are given in the Appendix, where gradients are also discussed.
    A brief summary of the approach is given in Algorithm \ref{alg:seq}. \resp{2-14}{The optimization with respect to $\tilde{\x}$ is non-convex; we pursue a multi-start approach by evaluating the acquisition function on a space filling design, then executing local searches via L-BFGS-B \citep{Morales2011} on the five points with the lowest acquisition function value.}

    \begin{algorithm}[ht]
	\caption{Pseudocode for the sequential learning approach}
	\begin{algorithmic}[1]
	\Require $n_0$, criteria $J(\cdot)$ (e.g., \texttt{Trace}, \texttt{Var1}, 
	\texttt{Var2})
	\State Construct and evaluate an initial design of experiments of size $n_0$ in 
	$\mathcal{X}$
	\State Build initial GP model
	\While{time/evaluation budget not exhausted}
	\State Find $\tilde{\x}^* \in \arg \max_{\mathbf{x} \in \mathcal{X}} J(\mathbf{x})$
	\State Evaluate $f(\tilde{\x}^*)$
	\State Update the GP model based on new data
	\EndWhile
	\end{algorithmic}
	\label{alg:seq}
	\end{algorithm}

    We now turn to application examples to demonstrate the benefits of
    sequential design for active subspace learning.
	

	\section{Numerical Experiments}
	\label{sec:experiments}
	
    We illustrate the advantage of judicious design-point selection by
    comparing sequential GP subspace estimation with random GP estimation and
    Monte Carlo estimation, as well as approaches using a global linear model
    (OLS) or local linear models (LL). The task at hand being estimation of a
    subspace, we measure error in a model's prediction as the sine of the
    first principal angle between the true active subspace and the active
    subspace estimated at each potentially expensive function evaluation. The
    sine of the first principal angle between two subspaces $\mathcal{U}$ and
    $\hat{\mathcal{U}}$, denoted as $\angle_1 {\mathcal{U}
    \hat{\mathcal{U}}}$, may be computed as the spectral norm of a simple
    expression of two unitary matrices, the first, $\mathbf{U}$, with range
    equal to $\mathcal{U}$ and the second, $\hat{\mathbf{U}}^\perp$, with
   
    kernel equal to $\hat{\mathcal{U}}$ \citep[Lemma 3.1]{Jiguang87}:
	\begin{equation}
	\label{eq:subspacemetric}
			 \sin\left(\angle_1 {\mathcal{U} \hat{\mathcal{U}}} \right) = 
\left\|\mathbf{U}^\top\hat{\mathbf{U}}^\perp\right\|_2.
	\end{equation}
	This ``subspace distance'' gives us some notion of the largest possible angle to be
	made between subspaces.
	
	All of the tested methods begin with the same
 	Latin hypercube sample and choose their next points either randomly or
 	based on an acquisition function.

	
	\subsection{The Rank-1 Quadratic}
	
	We begin by learning a quadratic function defined by a rank-1 matrix
	$\mathbf{A}$; that is,

	\begin{equation}
	\label{eq:rank1quad}
		f(\mathbf{x}) = \mathbf{x}^\top \mathbf{A} \mathbf{x} = \left(\mathbf{a}_1^\top\mathbf{x}\right)^2, \mathbf{x} \in [0,1]^m
	\end{equation}
    The gradient of this function is $\nabla f(\mathbf{x}) = 2\mathbf{A}\mathbf{x} = 2
    \left(\mathbf{a}_1^\top\mathbf{x}\right) \mathbf{a}_1$, so gradients live entirely
    in the one-dimensional subspace given by the range of $\mathbf{A}$.
   
   
    During each trial, a vector $\mathbf{a}_1$ is generated with iid standard normal entries.
    The results of applying each algorithm to 50 randomly generated rank-1
    $\mathbf{A}$ matrices
   
    in dimensions $m=$ 2, 5, and 8 are given in Figure~\ref{fig:quad}. The
    initial design is 5 times the dimension of the space, and the number of
    points selected sequentially is 10 times the dimension, giving a total
    budget of 15 times the input dimension in terms of function evaluations.
    To simulate numerical error or a stochastic experiment, we rerun the
    experiments this time adding Gaussian noise to the objective with standard deviation 
    $\tau = 5e^{-5}$. In the
    noiseless case, as soon as we see one gradient, we have the active
    subspace, so finite differencing (simple forward differencing with a step size of $10^{-4}$) together with the Monte Carlo estimator
    gets the right answer after $m+1$ evaluations. With slight numerical noise, however, 
    the accuracy obtained by finite differencing degrades, and the
    advantage of intelligent selection of evaluation points is clear.

	
	\begin{figure}[htpb]
		\centering
		\newcommand{0.39}{0.39}
		\includegraphics[scale=0.39]{images/quad/quad_2_d.pdf}
		\includegraphics[scale=0.39]{images/quad/quad_2_d_noisy.pdf}
		
		\includegraphics[scale=0.39]{images/quad/quad_5_d.pdf}
		\includegraphics[scale=0.39]{images/quad/quad_5_d_noisy.pdf}
		
		\includegraphics[scale=0.39]{images/quad/quad_8_d.pdf}
		\includegraphics[scale=0.39]{images/quad/quad_8_d_noisy.pdf}
		
		\caption{\small The quadratic problem \eqref{eq:rank1quad} in
		dimensions 2, 5, and 8. The $x$-axis represents computational effort
		in function evaluations while the $y$-axis represents the log of the mean (based on 50 trials) of the subspace error.
	 	The standard Monte Carlo estimator has zero error
		after a single gradient evaluation modulo minimal numerical noise in
		the deterministic case, so is not pictured in column 1, but would be at
		$-\infty$ in theory and was found to be around $10^{-10}$ in practice.}
		\label{fig:quad}
	\end{figure}


	The difference between the selected three definitions of variance of
	$\C^{(n+1)}$ can be observed in relation to the random infill, that is,
	randomly choosing the next point. \texttt{Var1} and \texttt{Var2} give
	similar results; \texttt{Trace} is in some cases barely better than random
	sampling. This indicates that to pinpoint the active subspace, it is insufficient to
	focus on the uncertainty associated with the mean (or sum) of the
	eigenvalues. Recall that each method begins with the exact same design, and thus that 
	discrepancies in accuracy at the beginning of the simulation are due to 
	how the methods differ in exploiting the same available data. Since the 
	classical MC estimator cannot avail itself of random function evaluations,
	it is given that many additional function evaluations to conduct finite differencing
	prior to the trials that appear in the figures.
	
	\subsection{The Wing Weight Function}
	
    The wing weight function, introduced by \citet{Forrester2011}, is a
    function of 10 inputs on a rectangular domain giving the weight of a light
    aircraft wing. Although it is a simple algebraic expression, it is of
    interest as a physically motivated test problem. Previously examined by
    the tutorial for the \texttt{active\_subspaces} \texttt{Python}
    package\footnote{\url{https://github.com/paulcon/active_subspaces/blob/master/tutorials/basic.ipynb}},
    a prominent one-dimensional active subspace was discovered \resp{1-1}{(see Figure \ref{fig:eigvals})}. Here, we
    explore this function beginning with a design of 20 points and choose an
    additional 80 points via a sequential design policy or randomly. Since an
    analytic representation of the active subspace is unavailable, we measure
    the ``true" active subspace by using finite differencing with the
    classical estimation technique with 10,000 function evaluations.
    Figure~\ref{fig:wing_weight} illustrates the advantage of design-point
    selection via our acquisition criteria for smaller function evaluation
    budgets. The variances that account for non-axis-aligned terms of
    $\C^{(n+1)}$ (i.e., in contrast to the \texttt{Trace} acquisition    
   
   
    function) perform best.
	
	\begin{figure}[ht]
		\centering
		\includegraphics[scale=0.45]{images/wing_weight_sim.pdf}
		\includegraphics[scale=0.45,trim={0 0 3.8cm 0},clip]{images/wing_weight_crits.pdf}
		\caption{Results on a simulation on the wing weight function. Left:
		Function evaluations are represented with the $x$-axis, while subspace
		error at each iteration is shown on the $y$-axis. The advantage of
		sequential design for smaller function evaluation budgets is clear.
		Right: The value maximizing the
		acquisition function is given for sequential strategies at each
		iteration.}
		\label{fig:wing_weight}
	\end{figure}

	\begin{figure}[ht]
		\centering
		\includegraphics[scale=0.45]{images/ww_eigs.pdf}
		\includegraphics[scale=0.45]{images/covid_eigs.pdf}
		\caption{\resp{}{The log of the eigenvalues of $\mathbf{C}$ matrices estimated with many data. Gaps indicate potential cut-offs for the active subspace. The Wing Weight active subspace matrix was estimated via Finite Differencing on 10,000 random points. The Covid-19 active subspace matrix was estimated via our GP methodology with 2,000 points, as the simulator's low accuracy precludes finite differencing. We chose active subspace dimensions of 1 and 3, respectively, as referenced in the experiments.}}
		\label{fig:eigvals}
	\end{figure}

	\resp{1-2}{
	\subsection{The Minnesota Covid-19 Simulator}
	
	The Minnesota Covid-19 Model \citep{mncovid} was developed to help guide public policy in response to the novel coronavirus, and calculates trajectories of various quantities of interest, such as cumulative death rates and ICU utilization over one year. However, the model depends on various epidemiological quantities, such as rate of asymptomatic infection, which are presently poorly understood. In this section, we treat this simulator as a black box and study the effect of varying certain parameters
	\footnote{Namely, the number of initial cases, the contact reduction of shelter in place and social distancing countermeasures, the proportion by which contact is reduced in persons aged 60 or more, the contact reduction due to personal behavior changes, the proportion of asymptomatic cases, the probability of a person over 80 being hospitalized, and the probability of a person over 80 dying at home}
	 on the total covid-19 deaths. The cumulative number of deaths, though calculated as a continuous variable, is not sufficiently precise so as to allow finite differencing. Thus, in our simulations, accuracy is computed as compared to the active subspace given by a Gaussian process fit on 2,000 randomly sampled data (Figure \ref{fig:eigvals} reveals that a three dimensional active subspace is reasonable). Figure \ref{fig:covid} illustrates the performance of our three sequential design criteria as compared to a random baseline. Again the variance-based acquisition functions perform best, while the \texttt{Trace} function is in this case only marginally better than random.
	
		\begin{figure}[ht]
		\centering
		\includegraphics[scale=0.45]{images/covid.pdf}
		\includegraphics[scale=0.45]{images/covid_boxplot.pdf}
		\caption{\resp{}{Comparison of GP based sequential design criteria on the Minnesota Covid-19 Simulator. Ground truth is taken to be the active subspace given by a GP with 2,000 randomly sampled input-output pairs. Left: Subspace error versus function evaluations, averaged over 40 runs. Right: Boxplots showing the final subspace error for each acquisition function over each of the 40 runs.}}
		\label{fig:covid}
	\end{figure}
}
	\section{Discussion and Perspectives}
	\label{sec:discussion}
	In this article, we have shown that the surrogate-assisted estimation of
	active subspaces, a technique previously executed via Monte Carlo finite
	differencing, has a closed-form solution. This result was leveraged to
	create acquisition functions enabling evaluation-efficient estimation of
	such subspaces via sequential design, potentially saving (expensive)
	black-box evaluations. Although in this article we initialized our

	search with a standard Latin hypercube sample, choosing an initial design
	related to our acquisition criteria may be preferable, an issue we leave
	for future work. Active subspace methods being a field of active research,
	we leave the estimation via Gaussian processes of recent refinements, such
	as that proposed by \cite{Lee2019}, for perspective. Future work could
	focus on problems specifically related to optimization, perhaps by
	considering a different measure to define the expected outer product of
	the gradient or by explicitly incorporating only those subspaces that play
	an important role in optimization, perhaps phrased in the framework of
	goal-oriented sensitivity analysis \citep{Fort2016}. Although not explored
	in this work, it is in principle straightforward to choose which points in
	the design space to explore in batches rather than one at a time. Such a
	problem, however, presents computational difficulty in optimization of the
	acquisition function. One may have to develop approximate methods in order
	for this to be carried out in practice.

	\section*{Acknowledgments}
	This material was based upon work supported by the U.S.\ Department of
	Energy, Office of Science, Office of Advanced Scientific Computing
	Research, applied mathematics and SciDAC programs under Contract No.\
	DE-AC02-06CH11357. The authors would like to thank Robert B.\ Gramacy for
	thoughtful comments on early drafts. \resp{}{This article benefited greatly from feedback provided by two anonymous referees.}
    
	