Deep active subspaces - a scalable method for high-dimensional uncertainty propagation

TL;DR

This paper introduces a scalable deep learning-based method for identifying active subspaces in high-dimensional uncertainty quantification problems, enabling efficient surrogate modeling without gradient information.

Contribution

It presents a novel approach combining reparameterization and deep neural networks to recover active subspaces without requiring gradient data, improving scalability.

Findings

Outperforms classical active subspaces in predictive accuracy

Effective in high-dimensional, data-scarce settings

Validated on synthetic and real-world datasets

Abstract

A problem of considerable importance within the field of uncertainty quantification (UQ) is the development of efficient methods for the construction of accurate surrogate models. Such efforts are particularly important to applications constrained by high-dimensional uncertain parameter spaces. The difficulty of accurate surrogate modeling in such systems, is further compounded by data scarcity brought about by the large cost of forward model evaluations. Traditional response surface techniques, such as Gaussian process regression (or Kriging) and polynomial chaos are difficult to scale to high dimensions. To make surrogate modeling tractable in expensive high-dimensional systems, one must resort to dimensionality reduction of the stochastic parameter space. A recent dimensionality reduction technique that has shown great promise is the method of `active subspaces'. The classical formulation of active subspaces, unfortunately, requires gradient information from the forward model - often impossible to obtain. In this work, we present a simple, scalable method for recovering active subspaces in high-dimensional stochastic systems, without gradient-information that relies on a reparameterization of the orthogonal active subspace projection matrix, and couple this formulation with deep neural networks. We demonstrate our approach on synthetic and real world datasets and show favorable predictive comparison to classical active subspaces.

Tables

Table 1. Table 1: Root mean square error (RMSE) on test dataset predictions from classic AS and deep AS response surfaces.

Approach	$\ell =1$	$\ell =0.01$
Classic AS	0.04378	0.17276
Deep AS	0.09241	0.23626

Equations

$$p(f)=\int \delta (f-f(\mathbf{\xi }))p(\mathbf{\xi })\mathrm{d}\mathbf{\xi },$$

$${\mu }_f=\int f(\mathbf{\xi })p(\mathbf{\xi })\mathrm{d}\mathbf{\xi },$$

$${\sigma }_f^2=\int (f(\mathbf{\xi })-{\mu }_f{)}^{2}p(\mathbf{\xi })\mathrm{d}\mathbf{\xi },$$

$$f(\mathbf{\xi })=g(\mathbf{\zeta })=g({\mathbf{W}}^T\mathbf{\xi }),$$

$${\mathcal{V}}_d\left({\mathbb{R}}^D\right)=\left\{\mathbf{A}\in {\mathbb{R}}^{D\times d}:{\mathbf{A}}^T\mathbf{A}={\mathbf{I}}_d\right\}.$$

$$\mathbf{C}=\int ({\nabla }_{\mathbf{\xi }}f(\mathbf{\xi }))({\nabla }_{\mathbf{\xi }}f(\mathbf{\xi }){)}^T\rho (\mathbf{\xi })\mathrm{d}\mathbf{\xi },$$

$$\mathbf{V}=\left[{\mathbf{V}}_1,{\mathbf{V}}_2\right],$$

$${\mathbf{C}}_S=\frac{1}{S}\sum _{i=1}^S{\mathbf{g}}^{(i)}{\mathbf{g}}^{(i),T}.$$

$$k_{\mathrm{AS}}(\mathbf{x},{\mathbf{x}}^{\prime })=k_{\mathrm{base}}(\mathbf{z},{\mathbf{z}}^{\prime })=k_{\mathrm{base}}({\mathbf{W}}^T\mathbf{x},{\mathbf{W}}^T{\mathbf{x}}^{\prime }),$$

$${\mathbf{W}}^{\ast },{\mathcal{H}}^{\ast },{\sigma }_n^{\ast }=\underset{\mathbf{W},\mathcal{H}}{\mathrm{argmax}}\log p(\mathbf{y}|\mathbf{X},\mathbf{W},\mathcal{H},{\sigma }_n),$$

$$\mathbf{W}=h(\mathbf{Q}),$$

$${\mathbf{w}}_i={\mathbf{w}}_{i-1}-\left(\frac{{\mathbf{w}}_{i-1}^T{\mathbf{q}}_i}{{\mathbf{w}}_{i-1}^T{\mathbf{w}}_{i-1}}\right){\mathbf{w}}_{i-1},i=2,3,\cdots ,d,$$

$$g(\mathbf{z})=f_{L+1}\circ f_L\circ \cdots \circ f_1(\mathbf{z}),$$

$$\hat{f}(\mathbf{\xi };\mathbf{\theta })=g(h(\mathbf{Q}{)}^T\mathbf{\xi };\mathbf{\theta }),$$

$$g(\mathbf{\zeta })={\mathbf{W}}^T\mathbf{\xi }=\alpha +{\mathbf{\beta }}^T\mathbf{z}+{\mathbf{z}}^T\mathbf{\Gamma }\mathbf{z}.$$

$$\nabla f(\mathbf{\xi })=\left(\beta +2{\mathbf{\xi }}^T\mathbf{W}\mathbf{\Gamma }\right){\mathbf{W}}^T.$$

$$\nabla \cdot \left(a(\mathbf{s})\nabla u(\mathbf{s})\right)=1,\mathbf{s}\in \Omega =[0,1{]}^2,$$

$$u(\mathbf{s})$$

$$\nabla u(\mathbf{s})\cdot \mathbf{n}$$

$$k(\mathbf{s},{\mathbf{s}}^{\prime })=\exp \left(-\frac{{\sum }_{i=1}^2|{\mathbf{s}}_i-{\mathbf{s}}_i^{{}^{\prime }}|}{\ell }\right),$$

$$\log a(\mathbf{s})=\sum _{i=1}^{100}\sqrt{{\lambda }_i}{\phi }_i(\mathbf{s})x_i,$$

$$q(\mathbf{x})=\mathcal{F}[u(\mathbf{s};\mathbf{x})]=\frac{1}{|{\Gamma }_2|}\int u(\mathbf{s};\mathbf{x})\mathrm{d}\mathbf{s}.$$

Full text

{nomenclature}\entry

$\mathbf{W}$Active subspace projection matrix \entry$g$Link function \entry$\boldsymbol{\xi}$Stochastic parameters \entryASActive subspace

1 Introduction

Inspite of the advent of the exascale era of computing and the rapid increase in the availability of computational resources [1], the sophistication of computer codes that simulate physical systems have also risen exponentially, either due to the incorporation of more realistic physics or higher-order numerical algorithms. A further consequence of the increasing sophistication of modern day computational solvers is the introduction of more parameters into the model to accurately describe boundary/initial conditions, material properties, constitutive laws etc. It is often the case that many (or all) of these parameters and unknown exactly. This brings up several questions for the computational scientist such as - i. how must one go about making robust predictions about the quantities of interest in a complex simulation, ii. how can one assess the impact of input parameter uncertainty on the model outputs, iii. how can one calibrate the model from experimental data, and so on. Answering such questions lie at the heart of uncertainty quantification (UQ) [2, 3].

The most common task in UQ is what is known as the forward UQ or uncertainty propagation (UP) problem [4, 5, 6]. A complete description of the UP problem can be found in Sec. 2.1. UP is the task of estimating the statistical properties of the model outputs given a formal description of the uncertainty in the model input parameters. Monte Carlo (MC) methods [7] are the most straightforward techniques for solving the UP problem and have, indeed, long been the workhorse of UQ [8, 9, 10]. A remarkable property of MC is that the variance of the standard MC estimate converges independent of the dimensionality of the stochastic parameters [7]. However, MC methods require a very large number of samples (hundreds of thousands) to show convergence in statistics [7]. This makes MC infeasible for state-of-the-art numerical simulators which have a large computational burden associated to each individual run.

To overcome the fact that one does not have the computational budget for hundreds of thousands of runs of a numerical simulator, one resorts to the surrogate approach to UP [11]. The idea is simple - perform a limited number of forward model evaluations, collect the resulting data, and construct a cheap-to-evaluate, yet accurate, map between the input uncertain parameters and the model outputs. This map serves as approximation to the true solver, and is referred to as a surrogate or a response surface. Since the response surface can be repeatedly evaluated very quickly, it is now easy to couple the surrogate with a MC approach to estimate model output statistics. Surrogate models can be either intrusive (i.e. requiring modification of the simulator for the analogous deterministic problem) or non-intrusive (i.e. where the application of the surrogate model is an ‘outer-loop’ process and the numerical solver can be treated as a black-box). Naturally, given the sophistication of state-of-the-art numerical models, intrusive methods such as polynomial chaos [12, 13] have fallen out of favor in recent times. This has coincided with the rise in popularity of non-intrusive methods such as Gaussian process regression (GPR) (or Kriging) [14, 15].

The surrogate approach to UQ has seen tremendous success in a broad range of applications [16, 17, 5]. However, state-of-the-art surrogate modeling techniques become exponentially difficult to scale as the number of stochastic parameters increases [18, 19]. This is due to the phenomenon known, universally, as the curse of dimensionality (CoD). Originally coined by Bellman in the context of dynamic programming [20], the CoD refers to the phenomenon where the volume of the input space that one must explore to gather data sufficient for constructing an accurate response surface, rises exponentially with a linear increase in the input dimensionality. To address the CoD, one needs to perform dimensionality reduction of the stochastic parameter space. The easiest approach to parameter space reduction is through a ranking of the importance of individual input dimensions. Methods that adopt this approach include sensitivity analysis [21] and automatic relevance determination (ARD) [22]. Unfortunately, such variable reduction techniques are most effective when the input variables have some degree of correlation. In generic UP problems, the stochastic parameters are frequently uncorrelated. For instance, the common scenario of functional uncertainties (such as random permeability in porous media flow) comprise of high (or infinite) dimensional stochastic parameter spaces with statistical independence between input dimensions. The most popular approach to dimensionality reduction within the UQ community is the truncated Karhunen-Loéve (KL) expansion [23]. The idea behind the truncated KL expansion is that a spectral decomposition of the uncertain parameters can be used to express them as a linear combination of an infinite number of iid (independent and identically distributed) random variables, following which the infinite series can be truncated by picking basis functions corresponding to the highest ranked eigenvalues by magnitude. This is the functional analogue of the well-known principal component analysis (PCA) [24] used extensively in the machine learning (ML) community. Although extensively used, the truncated KL expansion (or PCA) only considers information contained in the input data resulting in an overestimation of the intrinsic dimensionality of the stochastic parameter space.

Generally speaking, an effective technique for dimensionality reduction needs to exploit intrinsic structure within the underlying map being approximated. One such technique that has been recently popularized is the method of active subspaces (AS), introduced in [19]. An AS is a low-dimensional linear manifold embedded within the true high-dimensional parameter space, which maximally captures the variation of the underlying map. AS has been successfully applied to numerous applications since it’s introduction [25, 26, 27, 28]. However, a key drawback of the original AS framework is it’s reliance on gradient information about the model outputs, which are often difficult (if not impossible) to obtain. To overcome the gradient requirement for AS recovery, [18] proposed a framework which subsumes the AS projection matrix into the covariance kernel in GP regression and attempt to learn it from available data.

In this work, we propose a simple solution for recovering the AS and constructing surrogate models that does not rely on non-Euclidean algorithms for surrogate model optimization. Specifically, we express the AS projection matrix as the Gram-Schmidt orthogonalization of an unconstrained matrix. We rely on the fact that the Gram-Schmidt procedure is fully differentiable (and thus perfectly amenable to gradient-based learning through backpropagation). Furthermore, our method is agnostic to the specific class of function approximation. This is contrasted with the previous gradient-free AS framework proposed in [18]. Lastly, we couple this idea with deep neural networks (DNNs) and demonstrate true AS recovery. An added benefit of the proposed approach is that one can use it to improve DNN identifiability (even though it is not our primary concern).

This manuscript is organized as follows. We begin with a formal description of the UP problem and the surrogate approach to solving it in Sec. 2.1. We then, in Sec. 2, review the classical approach to AS recovery (see Sec. 2.2.1) and the gradient-free GPR based AS approach proposed in [18] (see Sec. 2.2.2). This is followed up with the description of our approach to AS, in Sec. 2.3. Finally, we demonstrate the proposed methodology on challenging high-dimensional surrogate modeling problems in Sec. 3 and demonstrate that we recover the true AS through comparisons with the classical approach.

2 METHODOLOGY

2.1 The surrogate approach to uncertainty propagation

Consider a physical system modeled with a (potentially complex, coupled) system of partial differential equations. The PDE(s) is solved numerically using a black-box computer code, which we denote as $f$. $f$ may be thought of as a multivariate function which accepts a vector of inputs $\boldsymbol{\xi}\in\boldsymbol{\Xi}\subset\mathbb{R}^{D}$ and produces a scalar quantity of interest (QoI) $f(\boldsymbol{\xi})\in\mathcal{Y}\subset\mathbb{R}$. Information about $f$ may be obtained through querying the solver at suitable input design locations $\boldsymbol{\xi}$. We allow for the possibility that our measurement from the computer code may be noisy, i.e., $y=f(\boldsymbol{\xi})+\epsilon$, where $\epsilon$ is a random variable (the measurement noise might arise as a consequence of quasi-random stochasticity or chaotic behavior). Given this setup, the uncertainty propagation (UP) task is summarized as follows. Given a formal description of the uncertainty in the input parameters, $\boldsymbol{\xi}\sim p(\boldsymbol{\xi})$, we would like to estimate the statistical properties of the QoI. This includes, the probability density,

[TABLE]

and measures of central tendency such as the mean:

[TABLE]

and variance:

[TABLE]

where, $\delta(\cdot)$ in Eqn. (1) refers to the Dirac $\delta$-function.

As already discussed in the introduction, the standard MC method is infeasible when there is a large computational cost associated with querying $f$ and one must resort to the surrogate approach - replacing the true simulator $f$ with an accurate, cheap-to-evaluate approximation, $\hat{f}$. To do this, one queries $f$ at a set of $N$ carefully selected design locations $\mathbf{X}=(\boldsymbol{\xi}^{i})_{i=1}^{N}$, resulting in a corresponding set of measurements, $\mathbf{y}=(y^{i})_{i=1}^{N}$. We refer to the observed data, collectively, as $\mathcal{D}=\{\mathbf{X},\mathbf{y}\}$. Although the task of careful selection of the input design locations are a subject of a great deal of research, an in-depth discussion of this topic is beyond the scope of the present work. Here we simply assume that we are given $\mathcal{D}$.

2.2 Active subspaces

The fact that we are working in a high-dimensional regime $D(\gg 1)$ makes the task of constructing an accurate surrogate model with limited data $\left(N\approx\mathcal{O}(D)\right)$ practically infeasible because of the curse of dimensionality. To circumvent this, one seeks to exploit low-rank structure within the true response $f$ and methods for doing so are broadly categorized as ‘dimensionality reduction’ techniques. In this work, we focus on the case where the response admits the following structure:

[TABLE]

where, $\mathbf{W}\in\mathbb{R}^{D\times d}$ is a tall-and-skinny matrix of orthogonal columns which projects the high-dimensional input $\boldsymbol{\xi}$ to $\boldsymbol{\zeta}$ lying in a $d$-dimensional subspace such that $d\ll D$. In particular, $\mathbf{W}$ is constrained to be an element of the set:

[TABLE]

$\mathcal{V}_{d}\left(\mathbb{R}^{D}\right)$ is known as the Stiefel manifold with $\mathbf{I}_{d}$ being the identity matrix in $\mathbb{R}^{d\times d}$ and $g:\mathbb{R}^{d}\rightarrow\mathbb{R}$ is known as the link function. The structure posited in Eqn. (4) takes on physical meaning where the columns of $\mathbf{W}$ correspond to directions of the input space most sensitive to variation in the response $f$. The dimensionality reduction induced by the introduction of this structure, significantly simplifies the task of learning an accurate surrogate model.

2.2.1 Classical approach to active subspaces

The classical approach to recovering the active subspace, introduced in [19], proceeds as follows. Let the gradient of the QoI w.r.t. the input be denoted as $\nabla_{\boldsymbol{\xi}}f=\left(\frac{\partial f}{\partial\boldsymbol{\xi}_{1}},\frac{\partial f}{\partial\boldsymbol{\xi}_{2}},\cdots,\frac{\partial f}{\partial\boldsymbol{\xi}_{D}}\right)\in\mathbb{R}^{D}$. Given a probability distribution $\rho$ endowed upon the input space, we define the symmetric positive semi-definite matrix,

[TABLE]

which admits the spectral decomposition $\mathbf{C}=\mathbf{V}\boldsymbol{\Lambda}\mathbf{V}^{T}$, where $\boldsymbol{\Lambda}$ is a diagonal matrix of eigenvalues ordered by magnitude. Separating the $d$ largest eigenvalues from the rest, we can write the matrix of eigenvectors, $\mathbf{V}$, as:

[TABLE]

where, $\mathbf{V}_{1}\in\mathbb{R}^{D\times d}$ is a matrix consisting of the eigenvectors corresponding to the $d$ largest eigenvalues and $\mathbf{V}_{2}\in\mathbb{R}^{D\times(D-d)}$ is composed of the remaining eigenvectors. The active subspace projection matrix, then, is simply, $\mathbf{W}=\mathbf{V}_{1}$. Since the integral in Eqn. (6) is intractable analytically (due to the black-box nature of $f$), one only has access to discrete samples of the gradient at input locations $\boldsymbol{\xi}$ sampled from the distribution $\rho$. Given a dataset of $S$ gradient evaluations, $\mathbf{g}^{(i)}=\nabla_{\boldsymbol{\xi}}f(\boldsymbol{\xi}^{(i)}),i=1,2,\cdots,S$, where the $\boldsymbol{\xi}^{(i)}$s are sampled iid from $\rho$, an approximation to the matrix $\mathbf{C}$ may be constructed as:

[TABLE]

One may think of the approximation $\mathbf{C}_{S}$ as an empirical covariance matrix. After recovering the projection matrix $\mathbf{W}$ through the above procedure, one can obtain projected inputs, using $\mathbf{z}=\mathbf{W}^{T}\mathbf{x}$ and using a suitable technique such as Kriging to learn the link function $g(\cdot)$.

2.2.2 Gradient-free approach to active subspaces

As discussed in Sec. 2.2.1, the classic approach to AS recovery requires the evaluation of an empirical covariance matrix from samples of the gradient $\nabla_{\boldsymbol{\xi}}f$. Obtaining gradient samples in challenging in practice. In some cases (such as simple dynamical systems), one might have access to an adjoint solver which can compute the gradients of the QoI wrt the input parameters [29]. In other cases, the gradients can be approximated through finite differences (FD). Note that a single first-order FD gradient evaluation requires 2 expensive forward model runs. Lastly, one might even approximate gradients through approximate global models for the data [30]. In general, the black-box nature of the response as well the associated cost of FD gradients means that one simply does not have access to $\nabla_{\boldsymbol{\xi}}f$ and therefore cannot construct $\mathbf{W}$ through the classical approach. To alleviate this limitation, [18] introduced a methodology for constructing surrogate models without requiring gradient information. The gradient-free approach relies on two key ideas:

In GPR, prior knowledge about the underlying function can be encoded in a principled manner through the mean and the covariance functions of the GP. Thus, a new covariance kernel may be defined where the AS projection matrix $\mathbf{W}$ is simply a hyperparameter and learned through available data, $\mathcal{D}$. Formally, the prior knowledge about the active subspace structure described in Eqn. (4) is expressed through a GP kernel which takes on the form:

[TABLE]

where, $k_{\mathrm{base}}(\cdot,\cdot)$ is any standard kernel (such as the Matern or Radial basis function (RBF) kernels) which expresses prior knowledge about the regularity properties of the link function $g(\cdot)$. Once the active subspace kernel has been defined, inference in GPR proceeds through the maximization of the log marginal likelihood of the data wrt the kernel hyperparameters i.e.,

[TABLE]

where, $\mathcal{H}$ is the set of all hyperparameters of the base kernel $k_{\mathrm{base}}$, and $\sigma_{n}$ is the standard deviation of the likelihood noise. 2. 2.

While it is easy to enforce positivity constraints on the hyperparameters $(\mathcal{H},\sigma_{n})=\boldsymbol{\phi}$, the optimization task in Eqn. (10) is made challenging because of the fact that it is non-trivial to enforce the orthogonality constraints on the projection matrix $\mathbf{W}$. In order to do so, the complete methodology of [18] relies on a coordinate-ascent scheme to iteratively optimize over the variables $\boldsymbol{\phi}$ while keeping $\mathbf{W}$ constant and vice versa. The optimization steps over $\boldsymbol{\phi}$ proceed via standard second-order techniques for unconstrained optimization, such as the L-BFGS method [31]. The optimization steps over the projection matrix $\mathbf{W}$ utilize an adapted version of gradient-ascent on the Stiefel manifold described in [32].

2.3 Deep active subspaces

The methodology introduced by [18] lifts the gradient requirement of the classical approach to AS recovery by subsuming the AS projection matrix into the covariance kernel of a GP. While the methodology is sound and experimentally shown to recover the true AS, it suffers from two major drawbacks -

It is not agnostic to the choice of the surrogate model. Note that the gradient-free method described in Sec. 2.2.2 necessitates a GP surrogate by construction. Inspite of the elegance of GPR, arising out of the principled framework it offers for incorporating prior knowledge, quantifying epistemic uncertainty and performing model selection, it’s standard formulation scales poorly due to the $\mathcal{O}(N^{3})$ inversion of the (potentially dense) covariance matrix required at each optimization step. While sparse GPR [33, 34] partially alleviates this poor scaling through the introduction of $M(\ll N)$ inducing variables or ‘pseudo-inputs’, the task of selecting or optimizing for the inducing input locations is non-trivial. 2. 2.

The proposed solution for optimizing over the projection matrix $\mathbf{W}$, while respecting orthogonality constraints, is itself non-trivial, introduces $Dd$ additional hyperparameters into the covariance kernel, and is prone to getting trapped in local stationary points [18].

We propose, here, a much simpler approach that is:

Is agnostic to the choice of the link function approximator, 2. 2.

Is trivial to implement.

Specifically, we express $\mathbf{W}$ as:

[TABLE]

where, $\mathbf{Q}\in\mathbb{R}^{D\times d}$ lies on the standard Euclidean space, and $h:\mathbb{R}^{D\times d}\rightarrow\mathbb{R}^{D\times d}$ orthogonalizes the columns of $\mathbf{Q}$. Specifically, we chose $h$ to be the celebrated Gram-Schmidt (GS) orthonormalization procedure [35]. The GS process may be summarized as follows. Given an unconstrained matrix $\mathbf{Q}=[\mathbf{q}_{1},\mathbf{q}_{2},\cdots,\mathbf{q}_{d}]\in\mathbb{R}^{D\times d}$, where the $\mathbf{q}_{i}$s are the columns of $\mathbf{Q}$, we apply the transformation,

[TABLE]

with $\mathbf{w}_{1}=\mathbf{q}_{1}$. The projection matrix $\mathbf{W}$ is then assembled by normalizing the $\mathbf{w}_{i}$s, i.e., $\mathbf{W}=\left[\frac{\mathbf{w}_{1}}{\|\mathbf{w}_{1}\|_{2}},\frac{\mathbf{w}_{2}}{\|\mathbf{w}_{2}\|_{2}},\cdots,\frac{\mathbf{w}_{d}}{\|\mathbf{w}_{d}\|_{2}}\right]$.

Now one only needs to care about the the Euclidean matrix $\mathbf{Q}$, and optimize it to the available data. Noting that the transformation specified by Eqn. (12) is fully differentiable (as it composed entirely of differentiable mathematical operations), one may simply define a routine implementing the GS process using a backpropagation-capable library (such as TensorFlow or PyTorch) to obtain exact gradients of any QoI wrt $\mathbf{Q}$.

Since the projection matrix $\mathbf{W}$ has been reparameterized without any concern for the specific structure of the link function, $g$, we are free to pick any suitable class of function approximator for $g$. In this work, we define $g$ to be a deep neural networks (DNN) [36], a class of highly flexible nonlinear function approximators with satisfy universal approximation properties. Formally, a $L$-layered DNN representation for $g$ is defined as:

[TABLE]

where, $f_{i}(\mathbf{z}_{i-1})=h_{i}(\mathbf{W}_{i}^{T}\mathbf{z}_{i-1}+\mathbf{b}_{i})$, with $\mathbf{W}_{i}\in\mathbb{R}^{d_{i}\times d_{i-1}},\mathbf{b}_{i}\in\mathbb{R}^{d_{i}},\mathbf{z}_{i}=\mathbb{R}^{d_{i}},\mathbf{z}_{0}=\mathbf{z},\mathbf{z}_{L+1}=g(\mathbf{z})$, and $h_{i}(\cdot)$ is a suitable nonlinear function applied elementwise on it’s argument. $h_{L+1}$ is set to be the identity function (since we are dealing with unconstrained real-valued outputs) and the other $h_{i}$s are set as the hyperbolic tangent function, a standard choice in the literature. The matrices $\mathbf{W}_{i}$s and the vectors $\mathbf{b}_{i}$s are called the ‘weights’ and ‘biases’ of the DNN and here we denote all of them collectively as $\boldsymbol{\theta}=\{\mathbf{W}_{1},\mathbf{W}_{2},\cdots,\mathbf{W}_{L+1},\mathbf{b}_{1},\mathbf{b}_{2},\cdots,\mathbf{b}_{L+1}\}$. The full surrogate, is there expressed as:

[TABLE]

where the unknown parameters $(\boldsymbol{\theta},\mathbf{Q})$ can be optimized through standard gradient-descent techniques. In this work, we use the famous Adaptive Moments (ADAM) optimization method [37].

3 NUMERICAL EXAMPLES

3.1 Synthetic example with known active subspace

Let $f:\mathbb{R}^{D}\rightarrow\mathbb{R}$ such that $f(\boldsymbol{\xi})=g(\boldsymbol{\zeta})=g(\mathbf{W}^{T}\boldsymbol{\xi})$ where $\mathbf{W}\in\mathcal{V}_{d}\left(\mathbb{R}^{D}\right)$. Define $g:\mathbb{R}^{d}\rightarrow\mathbb{R}$ as a quadratic function in $\mathbb{R}^{d}$:

[TABLE]

The gradients of $f$ are given by

[TABLE]

For this pedagogical example, we set $D=20$ and test our approach on two cases with true AS dimensionality, $d=1$ and $2$. The data for inputs $\boldsymbol{\xi}$, $\alpha$, $\beta$ and $\boldsymbol{\Gamma}$ are generated by sampling standard Gaussians of appropriate shapes. The matrix $\mathbf{W}$ is generated by performing the QR decomposition on a similarly generated matrix in $\mathbb{R}^{D\times d}$. Random seed is fixed for reproducibility. The output data $\mathbf{y}$ used for training is standardized, i.e. it is scaled to have 0 mean and unit variance.

3.1.1 Case 1: 1 dimensional active subspace

We begin testing our proposed gradient-free approach on a synthetic function which an AS of dimensionality $d=1$. To train the AS DNN, we use $N=50$ input-output observations. Furthermore, the output data is corrupted with zero mean Gaussian noise of standard deviation $1\times 10^{-2}$. The link function is approximated with a 2 layer DNN of 50 units per hidden layer and $L_{2}$ regularization with a weight decay constant of $1\times 10^{-4}$ is used to prevent overfitting. The ADAM optimizer is set to perform $3\times 10^{4}$ iterations with a base learning rate of $1\times 10{-3}$ dropped by a factor of $10$ every $10^{4}$ iterations. In Fig. 1 we visually compare the true AS for this case and the AS discovered by the DNN. We note that they are very close, indicating that our approach has found the correct AS. Fig. 1 also shows a comparison of the AS DNN predicted response with the true response from a test dataset of $500$ observations. Qualitatively, the predictions match the observations very closely. Quantitatively, we achieve a root mean square error of $0.039689$ on the test dataset. Note that we pursue no further optimization of our DNN structure.

3.2 Case 2: 2 dimensional active subspace

We now test our proposed on a synthetic function with AS dimensionality, $d=2$. We use $N=100$ input-output observations for training and corrupt the output data with zero mean Gaussian noise of standard deviation $1\times 10^{-2}$. We retain all other experimental settings from Sec. 3.1.1. A comparison of the true AS and the predicted AS shown in Fig. 2 reveals that we recover the low-dimensional quadratic response upto arbitrary rotations of the coordinate system. We compare the predicted response of the AS DNN and true outputs from a test dataset of $500$ observations. Again, we obtain excellent qualitative agreement as seen in Fig. 2 and quantitatively, we obtain a root mean squared error of $0.028748$ between the predicted and true outputs.

3.3 Benchmark elliptic PDE example

Consider the following stochastic elliptic partial differential equation defined on the unit square in $\mathbb{R}^{2}$:

[TABLE]

with boundary conditions:

[TABLE]

where, $\Gamma_{u}$ is the top, bottom and left boundaries and $\Gamma_{n}$ denotes the right boundary of $\Omega$. The diffusion coefficient $a$ (or conductivity field) is a spatially-varying uncertain input, and it’s logarithm is modeled as a 0-mean Gaussian random field, i.e., $\log a(\mathbf{s})\sim GP(a|0,k(\mathbf{s},\mathbf{s}^{\prime})))$, where, the covariance function $k$ is defined as follows:

[TABLE]

with $\ell$ being the correlation length. This formalization of the uncertainty around $a(\mathbf{s})$ makes it a stochastic process - an infinite dimensional quantity. We use the truncated KL expansion to perform a preliminary dimensionality reduction by expressing the logarithm of the field as:

[TABLE]

where, the $\lambda_{i}$s and the $\varphi_{i}$s are the eigenvalues and eigenfunctions of the correlation function, numerically obtained using the Nyström approximation [38], and the $x_{i}$s are uncorrelated, standard normal random variables. Denote all the $x_{i}$s collectively as $\mathbf{x}=(x_{1},x_{2},\cdots,x_{100})\sim\mathcal{N}(\mathbf{x}|\mathbf{0},\mathbf{I}_{100})$. We are interested in the following scalar QoI:

[TABLE]

Given a realization of the random variable, $\mathbf{x}=(x_{1},x_{2},\cdots,x_{1}00)$, one can generate a realization of the QoI, $q$, whose statistics one wishes to estimate. We have, at our disposal, a dataset of $300$ realizations of the random variable $\mathbf{x}$ and the corresponding solution $q$. With this dataset, we construct a surrogate that maps $\mathbf{x}$ to $q$, i.e., $\hat{f}:\mathbb{R}^{300}\rightarrow\mathbb{R}$. We will consider two cases of the correlation length $\ell$ - a short correlation length of $\ell=0.01$ and a long correlation length of $\ell=1$ and attempt to recover as AS with $d=1$. We randomly shuffle and split the data into a training set of $250$ samples, and test on the remaining $50$ samples. The output data is standardized to have zero mean and unit variance for numerical stability. The dataset for this example, and the code to generate it, can be found here: https://github.com/paulcon/as-data-sets/tree/master/Elliptic_PDE. Once again, we set our approximation of the link function to be a DNN with 2 hidden layers and 50 units per layer. All other experimental settings from Sec. 3.1 are retained. Lastly, for this example, samples of the QoI gradients are available and we use this to compare our results with the results obtained from classic AS. For the case of the classic AS, we use GPR as the link function.

Fig. 3 shows results comparing the deep AS approach and the classic AS approach for the $\ell=1$ case. Fig. 4 shows the same for the case of $\ell=0.01$. We observe that there is good qualitative agreement between the AS recovered by gradient-free deep AS approach and the gradient-based classic approach. This serves to verify the fact that our approach does indeed recover the correct AS. Table 1 shows a comparison of the RMSE error in the prediction of the outputs from the test dataset, for both $\ell$ cases. We observe that inspite of the fact that we do not use the information about the gradients of QoI, our gradient-free deep AS approach is are able to achieve RMSE comparable to the classic AS. Once again, we emphasize that we do not pursue optimization of the modeling choices involved in the DNN approximation of the link function.

An interesting observation that emerges from the comparison of the classic and deep AS approaches for the short correlation length in Fig. 4 is that inspite of recovering a one-dimensional AS, the test data predictions from both approaches show a discrepancy from their true values. Since the QoI is generated from a deterministic computer code, we cannot explain this deviation as ‘noise’. Rather, this suggests that a linear dimensionality reduction is sub-optimal and one might wish to recover a nonlinear generalization of the active subspace, such as the one discussed in [39]. An investigation into this shortcoming is beyond the scope of the present work. Finally, one may note that both the classic and the deep AS approaches perform much better on the $\ell=1$ case, relative to the $\ell=0.01$ case. This is unsurprising, considering that it becomes much more difficult to capture the uncertainty of the diffusion coefficient as it’s lengthscale reduces.

4 Conclusion

In this work, we presented a novel methodology for recovering active subspaces (AS) and constructing surrogate models in applications with high-dimensional uncertain parameter spaces. Our approach rests on a reparameterization of the AS projection matrix using the Gram-Schmidt procedure. Noting the fact that the GS procedure is, in principle, a fully-differentiable operation, one can easily backpropagate through the GS process to obtain gradients necessary in standard optimization routines. This formulation liberates us from the GPR approach of past gradient-free AS methods, and allows us to couple AS recovery with deep neural networks (DNNs) - a nonlinear function approximator which can be scaled to high-dimensions/larger datasets much more easily. We demonstrated the proposed approach on benchmark problems in AS recovery, showing that we do indeed recover the correct AS.

This work represents a first-step toward scaling gradient-free recovery of AS - an important objective since many quantities of interest encoding physical laws exhibit AS or AS-like ridge structure [40]. Our long term interest is in the development of efficient, Bayesian surrogates that are capable of quantify epistemic uncertainty. The reparameterization of the projection matrix allows us to leverage standard Bayesian inference methods such as stochastic variational inference (SVI) [41] or Hamiltonian Monte Carlo (HMC) [42] to construct Bayesian AS surrogates without resorting to specialized Riemannian manifold versions of these techniques.