\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.}


\section{Introduction}
Automatic Speech Recognition (ASR) field has been a research area focused on generating text from human speech. Traditionally, ASR system requires integrating different stages of models: namely acoustic, language and pronunciation models. An acoustic model first takes raw audio input and develop statistical models via extracting acoustic features such as Mel Frequency Cepstral Coefficients. These features are then applied to an observation model of the Hidden Markov Model (HMM) for temporal variability to estimate text sequence of the audio stream input. HMM integrates in its process the Language Model (LM) typically constructed of word sequence probabilities given a phone sequence \cite{Goodman2001}.

Over the last decade, there has been a dramatic shift in the ASR community due to the explosion of deep learning. Although the idea of deep layered network has been around over many decades, the sudden embrace of the method by the ASR community can be attributed to availability of both computational power and large public datasets that would allow training of deep layered networks possible.  Effectiveness of deep networks has been such that traditionally handcrafted feature based methods, such as GMM, was mostly replaced by the new paradigm. The approach yielded a significant improvement in performance of ASR systems \cite{Dahl2012a,Mohamed2012,Mohamed2009}.  Application of deep network to other areas such as lexicon and language followed for improved performance \cite{Rao2015,Mikolov2011,Xu2018,Sundermeyer2012}. In these earlier implementations of deep learning into ASR, the structures still require integration of these components mentioned.

Recently, end-to-end models for ASR have been an active research focus. These end-to-end models reduce the complexity of the system by making a direct mapping of input speech frames to output texts without the need of intermediate model such as GMM, HMM, or LM \cite{Hannun2014,Graves2014,Zhang2017b,Miao2015}. These end-to-end models, such as sequence-to-sequence \cite{Zhang2017b,Chan2016a,Bahdanau2015} or Connectionist Temporal Classification (CTC) based ASR \cite{Graves2014,Miao2015,Graves2006}, have achieved good performance in comparison to the conventional ASR models. These models reduce the complexity of deep learning based ASR systems significantly in terms of their structures as well as the training of the networks, and also delivers good performance.

As in most of the ASR approaches, including the ones mentioned here so far, however, the intelligibility of an ASR system is heavily influenced by presence of noise. While building an effective ASR system is a challenge, building an ASR system that is robust to noise remains more difficult. For traditional hand-crafted feature based systems, MFCC and log Mel-filter bank features are used as input to the model because they are based on the human auditory system. However, these features have no strategies for preventing the performance degradation in noisy conditions. For handling noise using deep learning architectures, there have been a variety of methods proposed including denoising autoencoder, which uses a non-linear mapping of noisy speech to target clean speech features. Deep Denoising Autoencoder (DDAE) \cite{Lu2013a,Feng2014a}, Denoising Variational Autoencoder (VAE) \cite{Hsu2018}, and multi-tasking autoencoder (MTAE) \cite{Zhang2018b} are some of these autoencoder based approaches applied to alleviate the noise issue.

We believe that ASR performance is highly dependent on how well noise problem is handled upfront to prevent corruption of signals for downstream chain. Thus, we propose a novel noise robust ASR approach with the following components. Within the front-end of the ASR framework, this work focuses on the feature enhancement method by appending the pre-processing stage by implementing a denoising autoencoder with skip connections where the original input is fed into different stages in the network. In addition, a correlation distance measure (CDM) that has been proposed in \cite{Szekely2007} and had been used in clustering \cite{JyotiBora2014} with better demonstrated performance over other similarity distance measures, is implemented.  The CDM is used in the objective function to increase the dependency between the latent and enhanced with the target features. The performance of our approach is measured by WER that is obtained by training a Recurrent Neural Network (RNN)-CTC phoneme-based ASR model using the EESEN toolkit .  We used DDAE \citep{Lu2013a}, and MTAE \cite{Zhang2018b} as baseline methods for comparison purpose.

\begin{figure*}[hbt!]
	\centering
	\@ifnextchar[{\ts@includegraphics}{\ts@includegraphics[width=\checkGraphicsWidth,height=\checkGraphicsHeight,keepaspectratio]}}[width=18cm, height=7cm]{Proposed_model_structure_v1_0.jpg}

	\caption{Skip-connection denoising autoencoder block diagram structure }
	\label{fig:1}
\end{figure*}

\section{Related Work}

\subsection{Denoising Autoencoder (DAE)} A denoising autoencoder, a variant of basic autoencoder, consists of two parts: the encoder and the decoder. The DAE is trained to reconstruct clean target speech from corrupted input. The encoder maps the input $\tilde{x} $ to a hidden representation $z=f(W\tilde{x}+b) $, and the decoder reconstructs an enhanced version $\hat{x}=f(Wz+b) $, where f denotes affine transformation \cite{Vincent2008}. This is done by minimizing the loss, computed from the mean squared error ($\textit{MSE} $) between the clean target ($\textit{x} $) and the output enhanced version of the model ($\hat{x}$),  as in Equation~(\ref{eq:1}).

\begin{equation}
\label{eq:1}
Loss(x,\hat{x})=||x-\hat{x}||^2
\end{equation}

where $\widetilde x $ corresponds to the corrupted speech features, $\hat{x}$ is the enhanced speech features, and x is the ground truth speech features \cite{Lu2013a,Feng2014a,Lu2012}. In a typical implementation, the speech features input to the DAE is not in the same dimension of the clean target speech features. This is because adjacent frames are also added as input by concatenating them with the current frame with knowledge that there is often leakage of speech information among neighbouring frames. Assuming noise being IID, however, contribution of noise from neighbouring frames should be minimal. More extended version of DAE has been proposed such as deep denoising autoencoder \cite{Lu2013a} with 5 hidden layers and using RBM training, and Multitasking Autoencoder where denoising task and noise features estimation task are performed separately in the DAE with mutually shared units among the two tasks \cite{Zhang2018b}.

\subsection{Skip connection networks} The skip connection, also known as shortcut connection, has been employed in image recognition and achieved successful results as in highway networks \cite{Srivastava2015,He2016b}.  It is also deployed in various image applications ranging from classification to image restoration by using autoencoder \cite{Mao2016}. In addition, the skip connection has been also used in acoustic field with applications in denoising autoencoder for music \cite{Liu2018}, speech enhancement \cite{Tu2017}, and others.

Additionally, different types of skip connection in image processing were investigated \cite{He2016} where, the identity shortcut had given the best performance. This identity shortcut has no extra parameters to be learned so it reduces the complexity and ensures information flow. Based on that, our proposed model is using an identity skip connection in the proposed structure.

\section{Proposed Model}
The proposed method is composed of two main parts: a novel skip connection architecture and a new objective function shown to be effective in noise removal. Novelty of our approach is an implementation of skip connections to the denoising autoencoder structure with an addition of a correlation distance measure in the objective function to enhance contribution of information from the skip connection for improved speech feature reconstruction in ASR. Figure~\ref{fig:1} illustrates the proposed skip-connection structure.

\subsection{Skip Connection Denoising Autoencoder (SK-DDAE) } \paragraph{Encoder} \sloppy The encoder extracts features of an input to a latent space by finding the best representation of the input with minimal influence by noise.  Input $\tilde{X}$ is fed to the proposed structure consisting of 11 concatenated frames of log Mel-filter bank ($\tilde{X}  = [\tilde{x}_{t-5},\tilde{x}_{t-4}, \dotsb ,\tilde{x}_{t-1},\tilde{x}_t,\tilde{x}_{t+1}, \dotsb, \tilde{x}_{t+4}, \tilde{x}_{t+5}] $) normalized in the range of [0,1] per utterance. We implement an identity skip connection in the encoder to allow the target frame, ($\tilde{x_t}$), to be fed directly into the middle layer of the encoder. The information of the skip connection ($\tilde{x_t}$) is concatenating with the output from the previous layer ($y_{i_enc-1}$), to be fed to the next layer of the encoder (i\textsuperscript{th} layer), and the output can be obtained as in Equation~(\ref{eq:2}).
\begin{equation}
\label{eq:2}
y_{i_{enc}}=\sigma(W.[y_{i_{enc}-1)},\tilde{x}_t ]+b)
\end{equation}

where $W$, $b$ and $\sigma$ are the weight, bias and the sigmoid function, respectively. The skip connection in our proposed model is somewhat different than the residual learning that was introduced in \citep{He2016b,Veit2016}. While their skip connection mechanism adds the input frame information, we instead concatenate the input frame $\tilde{x_t}$ to the mid-layer. We believe this approach improves the network to find correlated information in the required frame and carry it to influence the output for better representation in the output. The encoder in the proposed model consists of three layers of MLP with 512, 256 and 128 units.

\paragraph{Decoder} Purpose of the decoder is to reconstruct the target frame features from the latent features obtained from the encoder. The same skip connection (shortcut) from the input is also used as part of the decoder with the input frame concatenated to a mid-layer in the decoder based on the same premise that additional correlated information from the input frame can be gleaned in this stage. The decoder layers are organized in a symmetric manner to the encoder layers with 128, 256 and 512 used units, and finally only 40 enhanced features of the desired frame ($\hat{x}$) is obtained. The obtained features, eventually used for our ASR model, resulted better performance when noise is present in the input.

\subsection{Objective function}
In this paper, we are proposing a new objective function by adding a penalty term to the reconstruction error. The penalty is based on the Correlation Distance matrix approach that has been proposed in \cite{Szekely2007}. The correlation distance matrix had been used for clustering purpose as it is a similarity measure that considers as a true measure of dependency and test the joint independence between two random vectors. The dependency measure would add a new constraint in the network training to force the network to extract more salient features relevant to ASR.

\paragraph{Correlation Distance}  Assuming two random samples $(X,Y)={(X_k,Y_k ): k=1,\dotsb
, n} $ from a joint distribution of X in ${\mathbb{R}}^{n}$ and Y in ${\mathbb{R}}^{q}$, we compute the correlation  distance measure($\mathcal{\Re}(X,Y)$) as in Equation~(\ref{eq:3}).

\begin{equation}
\label{eq:3}
\begin{aligned}
\mathcal{\Re}^2(X,Y)=
\begin{cases}
\frac{\mathcal{V}^2(X,Y)}{\sqrt[2]{\mathcal{V}^2(X)\mathcal{V}^2(Y)}}, 	&\text{$\mathcal{V}^2(X)\mathcal{V}^2(Y)>0$,}\\

0,		&\text{$\mathcal{V}^2(X)\mathcal{V}^2(Y)=0$}
\end{cases}
\end{aligned}
\end{equation}

where,  $V^{2}(X,Y) $is the covariance distance that is obtained by computing the double centered distance matrices A and B by defining:

\begin{equation}
\label{eq:4}
\begin{aligned}
a_{kl}=|X_k- X_l|_p ,	\bar{a}_{k.}=1/n\sum_{l=1}^{n} a_{kl} ,\\
	\bar{a}_{.l}=1/n\sum_{k=1}^{n}a_{kl} ,	\bar{a}_{..}=1/n^2 \sum_{k,l=1}^n a_{kl}
	\end{aligned}
\end{equation}
\begin{equation}
\label{eq:5}
A_{kl}=a_{kl}-\bar{a}_{k.}-\bar{a}_{.l}+\bar{a}_{..}
\end{equation}

Similarly,
\begin{equation}
\begin{aligned}
\label{eq:6}
b_{kl}=|Y_k- Y_l |_p ,	\bar{b}_{k.}=1/n \sum_{l=1}^nb_{kl} ,\\
	\bar{b}_{.l}=1/n\sum_{k=1}^nb_{kl} ,	\bar{b}_{..}=1/n^2\sum_{(k,l=1)}^{n}b_{kl}
\end{aligned}
\end{equation}

\begin{equation}
\label{eq:7}
B_{kl}=b_{kl}-\bar{b}_{k.}-\bar{b}_{.l}+\bar{b}_{..}
\end{equation}

Hence,

\begin{equation}
\label{eq:8}
\begin{aligned}
V^2 (X,Y)=1/n^2\sum_{(k,l=1)}^{n}{A_{kl} B_{kl}}
\end{aligned}
\end{equation}

The  $\mathcal{\Re}(X,Y) $ guarantees nonnegativity and it is in range of [0,1] in contrast to the Pearson correlation coefficient that has a range of [-1,1]. Additionally, when  $\mathcal{\Re}(X,Y) = 0 $, it indicates independence between these two random vectors. However, when $\mathcal{\Re}(X,Y) = 1 $, it indicates linear dependency \cite{Vincent2008}.

Therefore, the objective function that is used for training the SK-DDAE can be expressed as:

\begin{equation}
\begin{aligned}
\label{eq:9}
loss=\norm{x_t-\hat{x}}_2^2+\beta*((1-\mathcal{\Re}(z,x_t))+\\
(1-\mathcal{\Re}(\hat{x},x_t)))
\end{aligned}
\end{equation}

  where, $z$ and $\beta$ are the latent features that are obtained from the encoder part and penalty term parameter, respectively. The objective function will maximize the dependency between the latent features and the target, and would drive the network toward producing $z$ that has good representation of the target by minimizing noise information. Additionally, it also measures the dependency between the output and target, and try to increase the dependency in the training process. Furthermore, a second variation of the proposed model was investigated by using a nonlinear term in the objective function (Energy term) that we assume will enhance the system performance as:

\begin{multline}
\label{eq:10}
loss=\norm{x_t-\hat{x}}_2^2+\beta*(1-\mathcal{\Re}(z,x_t )+(1-\mathcal{\Re}(\hat{x},x_t ))\\
+\sigma*((1-\mathcal{\Re}(z,x_t))^2+(1-\mathcal{\Re}(\hat{x},x_t))^2
\end{multline}

In this paper, we will refer the model that uses only the mean squared error loss as SK-DAE, the linear penalty correlation distance term as CDSK-DAE, and the nonlinear penalty correlation distance term as CDESK-DAE.

In addition, we use $\beta$, $\sigma $ = 0.01\$ for simplicity, which we leave the process of finding the optimal penalty term parameters ($\beta$ and $\sigma$) for future work. Also, we found that the structure performs better with identity skip connection than using CNN skip connection which is consistent with the finding in \cite{He2016}.

\begin{table*}[!htbp]
\caption{{Results comparison by using word error rate (WER\%) per method in clean condition} }
\label{1}
\def1{1}
\ignorespaces
\centering
\begin{tabulary}{\linewidth}{LLLLLL}
\hlinewd{.8pt}\\[-12pt]
Log Mel-Filter Bank &  MTAE (WER\%) &  DDAE (WER\%) &  SK-DDAE (WER\%)  &  CDESK-DAE (WER\%)  &  CDSK-DAE (WER\%) \\
7.92 &  9.85 &  10.01 &  9.54  &  8.49  &   8.65 \\
\noalign{\vspace*{6pt}}\hline\noalign{\vspace*{2pt}}
\end{tabulary}\par
\end{table*}

\begin{table*}[!htbp]
\caption{{Results comparison by using word error rate percentage (WER\%) between the baselines and proposed model against different type of seen noise and levels of SNR} }
\label{2}
\def1{1}
\ignorespaces
\centering
\begin{tabulary}{\linewidth}{LLLLLLLL}
\hlinewd{.8pt}\\[-12pt]
SNR (dB)  &  Mel Filter Bank (WER\%)  &  MTAE (WER\%)  &  DDAE (WER\%)  &  SK-DDAE (WER\%)  &  CDESK-DAE (WER\%) &  CDSK-DAE (WER\%) &  \\
\multicolumn{8}{p{\dimexpr(\mcWidth{1}+\mcWidth{2}+\mcWidth{3}+\mcWidth{4}+\mcWidth{5}+\mcWidth{6}+\mcWidth{7}+\mcWidth{8})}}{Shopping center }\\
0 &  55.18 &  40.14 &  44.59 &  39.09 &  37.46 &  38.53 &  \\
5 &  35.58 &  20.96 &  28.5 &  25.5 &  25.59 &  27.7 &  \\
  10 &  20.27 &  23.27 &  19.74 &  18.16 &  18.68 &  19.32 &  \\
  20 &  10.81 &  15.15 &  13.59 &  12.62 &  12.51 &  12.92 &  \\
\multicolumn{8}{p{\dimexpr(\mcWidth{1}+\mcWidth{2}+\mcWidth{3}+\mcWidth{4}+\mcWidth{5}+\mcWidth{6}+\mcWidth{7}+\mcWidth{8})}}{Pub}\\
0 &  90.31 &  65.53 &  78.77 &  70.39 &  71.54 &  69.36 &  \\
5 &  60.13 &  48.56 &  51.02 &  46.73 &  46.59 &  45.74 &  \\
10 &  31.84 &  32.84 &  30.52 &  29.24 &  28.51 &  28.65 &  \\
20 &  11.84 &  17.15 &  15.13 &  14.71 &  14.67 &  13.96 &  \\
\multicolumn{8}{p{\dimexpr(\mcWidth{1}+\mcWidth{2}+\mcWidth{3}+\mcWidth{4}+\mcWidth{5}+\mcWidth{6}+\mcWidth{7}+\mcWidth{8})}}{Caf{\'e}}\\
0 &  74.87 &  56.09 &  62.98 &  55.18 &  54.03 &  55.38 &  \\
5 &  51.73 &  41.66 &  43.95 &  39.04 &  37.34 &  40.17 &  \\
10 &  30 &  32.27 &  29.72 &  27.27 &  27.04 &  29.22 &  \\
20 &  11.78 &  18.08 &  15.88 &  14.55 &  15.35 &  15.01 &  \\
\multicolumn{8}{p{\dimexpr(\mcWidth{1}+\mcWidth{2}+\mcWidth{3}+\mcWidth{4}+\mcWidth{5}+\mcWidth{6}+\mcWidth{7}+\mcWidth{8})}}{Schoolyard}\\
0 &  82.79 &  61.79 &  66.77 &  62.61 &  59.4 &  60.22 &  \\
5 &  56.23 &  43.15 &  42.71 &  39.06 &  36.75 &  38.88 &  \\
10 &  33.19 &  31.35 &  25.66 &  25.5 &  24.76 &  27.08 &  \\
20 &  11.93 &  17.17 &  13.72 &  13.91 &  13.91 &  13.2 &  \\
\noalign{\vspace*{6pt}}\hline\noalign{\vspace*{2pt}}
\end{tabulary}\par
\end{table*}

\section{Experiment}
The experiment is conducted in two stages; first one is about training the skip- connection denoising autoencoder to obtain the enhanced features. In the second stage, the enhanced features are used in the CTC ASR model based on the EESEN to evaluate the performance of the proposed method.

\subsection{Dataset}The skip-connection autoencoder is trained using the TIMIT dataset, which uses 3696 utterances of clean speech and 14784 augmented utterances with background noise.  The TIMIT dataset (3696 utterance) contains 3.14 hours of data and it is an Acoustic- Phonetic Continuous Speech Corpus, which has 61 phonemes in total. The TIMIT utterances were augmented with 4 types of noise (caf{\'e}, pub, schoolyard, shopping center) with four different SNR levels (0dB, 5dB, 10dB, 20dB). On the other hand, the EESEN ASR model was trained on the WSJ's si-84 data subset which was divided into 95\% for training and 5\% for validation. For evaluation, the WSJ test set (Nov 92) data is augmented with seven type of seen and unseen noise for the SK-DDAE (caf{\'e}, pub, schoolyard, shopping center, living room, nature creek, outside traffic road) obtained from the ETSI background noise database \cite{EuropeanTelecommunicationsStandardsInstitute2008} using ADDNOISE MATLAB library \cite{ITU-T2011a}. In the evaluation, background noises had been added with four different SNR level.s (0dB, 5dB, 10dB, 20dB). Both the WSJ and TIMIT data are obtained from Linguistic Data Consortium (LDC).

\subsection{Experimental setting}The log Mel-filter bank features were extracted with 512 FFT bins with 40 filterbanks to obtain 40 dimensional features. These features are normalized using CMVN and then the delta-delta is added to create a total of 120 features to train the acoustic model.  Furthermore, the denoising autoencoder baselines (DDAE and MTAE) and the proposed methods, are trained using normalized log Mel-filterbank features with batch size of 500. Both the DDAE and MTAE are trained with 50 epochs, whereas, the SK-DDAE is trained with 16 epochs. The proposed model is initialized using the Xavier initialization \cite{Jia2014} and uses AdamOptimizer with learning rate of 0.001. For Evaluation, normalized 40 dimensional features per frame was obtained for the WSJ dataset. A delta-delta features are added to the enhanced features to obtain 120 dimensions that are used for the acoustic model. Furthermore, the acoustic model is a RNN-CTC phoneme based model that have been developed using the EESEN toolkit and uses the same configuration as described in \cite{Miao2015}. The RNN is 4 bidirectional LSTM layers with 320 memory cells for both forward and backward sub-layers.

\begin{table*}[!htbp]
\caption{{Results comparison by using word error rate percentage (WER\%) between the deep learning baselines and proposed model against different type of unseen noise and levels of SNR } }
\label{3}
\def1{1}
\ignorespaces
\centering
\begin{tabulary}{\linewidth}{LLLLLL}
\hlinewd{.8pt}\\[-12pt]
SNR (dB) &
  MTAE (WER\%) &
  DDAE (WER\%) &
  SK-DDAE (WER\%) &
  CDESK-DAE (WER\%) &
  CDSK-DAE (WER\%)\\
\multicolumn{6}{p{\dimexpr(\mcWidth{1}+\mcWidth{2}+\mcWidth{3}+\mcWidth{4}+\mcWidth{5}+\mcWidth{6})}}{Livingroom }\\
0 &  60.3 &  53.41 &  50.77 &  50.54 &  49.78\\
5 &  42.94 &  36.43 &  35.87 &  34.95 &  35.25\\
10 &  30.11 &  24.03 &  24.26 &  23.87 &  25.31\\
20 &  15.97 &  14.21 &  13.37 &  13.29 &  13.63\\
\multicolumn{6}{p{\dimexpr(\mcWidth{1}+\mcWidth{2}+\mcWidth{3}+\mcWidth{4}+\mcWidth{5}+\mcWidth{6})}}{Natural Creek }\\
0 &  95.16 &  92.79 &  87.17 &  88.06 &  86.85\\
5 &  70.88 &  70.81 &  65.98 &  62.93 &  63.76\\
10 &  45.83 &  43.03 &  40.44 &  38.6 &  40.23\\
20 &  19.74 &  16.06 &  15.77 &  16.27 &  15.59\\
\multicolumn{6}{p{\dimexpr(\mcWidth{1}+\mcWidth{2}+\mcWidth{3}+\mcWidth{4}+\mcWidth{5}+\mcWidth{6})}}{Outside Traffic Road }\\
0 &  67.91 &  66.12 &  61.01 &  60.62 &  60.77\\
5 &  47.26 &  44.21 &  40.44 &  41.02 &  41.36\\
10 &  32.02 &  26.72 &  24.21 &  24.3 &  25.22\\
20 &  15.9 &  13.98 &  13.7 &  13.88 &  13.66\\
\noalign{\vspace*{6pt}}\hline\noalign{\vspace*{2pt}}
\end{tabulary}\par
\end{table*}

\begin{table*}[!htbp]
\caption{{Correlation distance measure among the noise and speech utterance in TIMIT training data} }
\label{4}
\def1{1}
\ignorespaces
\centering
\begin{tabulary}{\linewidth}{p{4.5em}p{4.5em}p{4.5em}p{4.5em}p{4.5em}p{4.5em}p{4.5em}p{4.5em}p{4.5em}}
\hlinewd{.8pt}\\[-12pt]
    & Shopping Center &  Caf{\'e} &  Pub &  Schoolyard  &
  Livingroom  &  Natural Creek &  Outside Traffic Road\\
Speech &  0.1473	& 0.1441	& 0.1781	& 0.1897	& 0.1419	& 0.1435  &	0.1699 \\
\noalign{\vspace*{6pt}}\hline\noalign{\vspace*{2pt}}
\end{tabulary}\par
\end{table*}

\begin{table}[!htbp]
\caption{{Correlation distance measure among the training and testing background noise} }
\label{5}
\def1{1}
\ignorespaces
\centering
\begin{tabulary}{\linewidth}{p{4.3em}p{4.3em}p{4.3em}p{4.3em}p{4.3em}
\hlinewd{.8pt}\\[-12pt]
Noise Type	& Livingroom	& Natural Creak	& Outside Traffic Road\\
Shopping Center & 0.2161 & 0.2249 & 0.2542\\
Caf{\'e} & 0.2042 & 0.1850 & 0.2664\\
Pub & 0.2271 & 0.2439 & 0.2576\\
Schoolyard & 0.2597 & 0.2515 & 0.3041\\
\noalign{\vspace*{6pt}}\hline\noalign{\vspace*{2pt}}
\end{tabulary}\par
\end{table}

\section{Results and Discussion}
The proposed model was tested in three different sets of values for $\beta$ and $\sigma$ ($\beta$=0, $\sigma$=0 for SK-DDAE, $\beta$=0.01 and $\sigma$=0 for  CDSK-DAE, and both $\beta$, $\sigma$=0.01 for CDESK-DAE). The result on Table~\ref*{1} illustrates the word error rate on the WSJ database in a clean environment. Both CDESK-DAE and CDSK-DAE performed better than the other deep learning methods (DDAE and MTAE) in identical conditions.

In both Tables~\ref*{2} and~\ref*{3} , the results illustrate the performance sensitivity of the baseline and the proposed model to various levels of SNR, as well the types of noise, both seen and unseen, in Word Error Rate (WER). Based on these results, the proposed models (CDSK-DAE and CDESK-DAE) outperformed the baseline models (MTAE and DDAE) under all noise types and under all SNR for seen and unseen noises except for very few cases. The performance improvement rates of CDESK-DAE ranging from 7.25\% to 10.14\% for seen noises and 13.91\% and 6.67\% for unseen noise were obtained when the performance of MTAE and DDAE were used as baselines. Similarly, CDSK-DDAE demonstrated  improvement rates ranging 5.27\% to 8.21\% for seen noise and 13.35\% to 6.06\% for unseen noise. Even SK-DDAE which is only trained with mean squared error showed improvement rates in comparison to MTAE and DDAE for seen noise ranging 5.59\% to 8.52\% and unseen noise ranging 13.06\% to 5.74\%, respectively. More significantly, the proposed models have demonstrated generalization capability by showing significant performance improvements of the WER\% in the cases of the unseen noises.

Figure \ref{fig:2} demonstrates the correlation distance matrix of the models (SK-DAE, CDMSK-DAE and CDMESK-DAE) during the training between the latent features with the target features. The addition correlation distance measure in the objective function increases the dependency or relation between the extracted latent features and the target features as shown in the CDMSK-DAE and CDMESK-DAE, in contrast to the SK-DDAE. This increase of dependency results in obtaining a better representation of our target by extracted relevant features in the latent space that helps in the reconstruction process. Similarly, Figure \ref{fig:3} shows that the correlation distance matrix of these three models during the training between the obtained enhanced features ($\hat{x}$) and the target features ($x$). Where the proposed objective function in both CDSK-DAE and CDESK-DAE gives slightly stable high dependency that is maintained through the training than the SK-DDAE in the first epochs.

\begin{figure}
	\centering
	\@ifnextchar[{\ts@includegraphics}{\ts@includegraphics[width=\checkGraphicsWidth,height=\checkGraphicsHeight,keepaspectratio]}}{CMD_latent_target_features_v1.png}
	\vspace*{-5mm}
	\caption{Distance correlation matrix during training between the latent features ($z$) and target clean features ($x$)}
	\label{fig:2}
\end{figure}

\begin{figure}
	\centering
	\@ifnextchar[{\ts@includegraphics}{\ts@includegraphics[width=\checkGraphicsWidth,height=\checkGraphicsHeight,keepaspectratio]}}{CMD_enhanced_target_features_v1.png}
	\vspace*{-5mm}
	\caption{Distance correlation matrix during training between the enhanced features ($\hat{x}$) and the target clean features ($x$)}
	\label{fig:3}
\end{figure}

Furthermore, Table \ref{4} illustrates the dependency among both the speech utterance (TIMIT train dataset) and the noise type that is used in this work using distance correlation measure. Pub, schoolyard and outside traffic road background noise have the highest correlation distance measure and that is reflected also from the high obtained WER\% from the ASR model. Although, natural creek background noise was not much correlated with the speech, it still gives the highest WER\%. That is because it was less correlated with the training background noise as illustrated in Table \ref{5}.

In summary, the results show that the proposed model with the non-linear correlation distance measure penalty (CDESK-DAE) has the best overall performance under seen and unseen noisy environment. Additionally, the obtained results demonstrate effectiveness of the objective function that uses the correlation distance measure. As assumed, the correlation distance measure seems to increase such dependency between the model and the target during the training, which resulted in better enhanced features that can be used for the ASR system. Although the newly introduced objective function clearly demonstrates its effectiveness, we believe that further research could prompt additional improvements for the denoising task in the ASR applications and beyond.

\section{Conclusions}
In this paper, we proposed a novel architecture of denoising autoencoder by developing a new objective function that uses a correlation distance measure to enhance the dependency between the model and the target features. The model was able to enhance the corrupted features with noise and showed an improvement of performance in both seen and unseen noisy environments. Although the models were effective, we believe further performance improvement can be achieved by optimizing the parameters $\beta$ and $\sigma$.  Additionally, we will investigate the objective function's effectiveness in other acoustic applications.

\section*{Acknowledgments}The work has been supported by the Korea Environmental Industry and Technology Institute's environmental policy-based public technology development project (2017000210001) that is funded by the Ministry of Environment and the contribution of David was supported by the US Army Research Laboratory.

\bibliographystyle{vancouver}