Extremum Sensitivity Analysis with Polynomial Monte Carlo Filtering

TL;DR

This paper enhances global sensitivity analysis by introducing efficient computation methods for sensitivity indices using polynomial ridge approximations and proposes a novel extremum Sobol' index based on Monte Carlo filtering, with empirical exploration of their relationship.

Contribution

It introduces a new approach to compute sensitivity indices efficiently with polynomial ridge approximations and proposes extremum Sobol' indices based on Monte Carlo filtering.

Findings

Sensitivity indices can be computed with fewer model evaluations using ridge approximations.

Extremum Sobol' indices based on Monte Carlo filtering are effective near output extrema.

Empirical analysis shows a relationship between MCF-based indices and skewness-based indices.

Abstract

Global sensitivity analysis is a powerful set of ideas and heuristics for understanding the importance and interplay between uncertain parameters in a computational model. Such a model is characterized by a set of input parameters and an output quantity of interest, where we typically assume that the inputs are independent and their marginal densities are known. If the output quantity is smooth, polynomial chaos can be used to extract Sobol' indices. In this paper, we build on these well-known ideas by examining two different aspects of this paradigm. First, we study whether sensitivity indices can be computed efficiently if one leverages a polynomial ridge approximation -- a polynomial fit over a subspace. Given the assumption of anisotropy in the dependence of a function, we show that sensitivity indices can be computed with a reduced number of model evaluations. Second, we discuss methods for evaluating sensitivities when constrained near output extrema. Methods based on the analysis of skewness are reviewed and a novel type of indices based on Monte Carlo filtering (MCF) -- extremum Sobol' indices -- is proposed. We combine these two ideas by showing that these indices can be computed efficiently with ridge approximations, and explore the relationship between MCF-based indices and skewness-based indices empirically.

Tables

Table 1. Table 1: Input variables to the borehole function. Distribution parameters to a Gaussian distributed variable refer to the mean and standard deviation respectively; for a uniformly distributed variable, they refer to the lower and upper limits of the support.

Input	Distribution	Distribution parameters	Description
$r_w$	Gaussian	0.10, 0.0161812	Radius of borehole (m)
$r$	Uniform	100, 50000	Radius of influence (m)
$T_u$	Uniform	63070, 115600	Transmissivity of upper aquifer (m2/yr)
$H_u$	Uniform	990, 1110	Potentiometric head of upper aquifer (m)
$T_l$	Uniform	63.1, 116	Transmissivity of lower aquifer (m2/yr)
$H_l$	Uniform	700, 820	Potentiometric head of lower aquifer (m)
$L$	Uniform	1120, 1680	Length of borehole (m)
$K_w$	Uniform	1500, 15000	Hydraulic conductivity of borehole (m/yr)

Table 2. Table 2: Input variables to the piston model.

Input	Range	Description
$M$	[30, 60]	Piston mass (kg)
$S$	[0.005, 0.020]	Piston surface area (m2)
$V$	[0.002, 0.010]	Initial gas volume (m3)
$k$	[1000, 5000]	Spring coefficient (N/m)
$P_0$	[90000, 110000]	Atmospheric pressure (N/m2)
$T_a$	[290, 296]	Ambient temperature (K)
$T_0$	[340, 360]	Filling gas temperature (K)

Equations

$$p=p\left({\mathbf{M}}^T\mathbf{x}\right)\text{where}\mathbf{M}\in {\mathbb{R}}^{d\times n}\text{and}\mathbf{x}\in {\mathbb{R}}^d.$$

$${\left({\mathbf{x}}_i,f_i\right)}_{i=1}^N,$$

$$\mathbf{f}=\left[\begin{array}[]{c}f_{1}\\ \vdots\\ f_{N}\end{array}\right]\approx\left[\begin{array}[]{c}p_{d}\left(\mathbf{x}_{1}\right)\\ \vdots\\ p_{d}\left(\mathbf{x}_{N}\right)\end{array}\right]=\bm{A}_{d}\bm{\alpha},$$

$$p_d(\mathbf{x})=\sum _{i=1}^r{\alpha }_i{\Psi }_i(\mathbf{x}),$$

$${\mathbf{A}}_d=\left[\begin{array}{cccc}{\Psi }_1({\mathbf{x}}_1)& {\Psi }_2({\mathbf{x}}_1)& \dots & {\Psi }_r({\mathbf{x}}_1)\\ {\Psi }_1({\mathbf{x}}_2)& {\Psi }_2({\mathbf{x}}_2)& \dots & {\Psi }_r({\mathbf{x}}_2)\\ ⋮& ⋮& \ddots & ⋮\\ {\Psi }_1({\mathbf{x}}_N)& {\Psi }_2({\mathbf{x}}_N)& \dots & {\Psi }_r({\mathbf{x}}_N)\end{array}\right].$$

$$\underset{\mathbf{\alpha }}{\mathrm{minimize}}$$

$$\begin{array}{rl}\underset{\mathbf{\alpha }}{\mathrm{minimize}}& {\Vert \mathbf{\alpha }\Vert }_1\\ \text{subject to}& {\Vert {\mathbf{A}}_d\mathbf{\alpha }-\mathbf{f}\Vert }_2\le ϵ,\end{array}$$

$$\mathbf{f}=\left[\begin{array}[]{c}f_{1}\\ \vdots\\ f_{N}\end{array}\right]\approx\left[\begin{array}[]{c}p_{n}\left(\bm{M}^{T}\mathbf{x}_{1}\right)\\ \vdots\\ p_{n}\left(\bm{M}^{T}\mathbf{x}_{N}\right)\end{array}\right]=\bm{A}_{n}\bm{\beta},$$

$${\mathbf{A}}_d\mathbf{\alpha }={\mathbf{A}}_n\mathbf{\beta },$$

$$\rho (\mathbf{x})=\prod _{j=1}^d{\rho }_j(x_j).$$

$${\sigma }_S=\frac{\text{Var}[\mathbb{E}[f(\mathbf{x})|{\mathbf{x}}_S]]}{\sigma },$$

$$f(\mathbf{x})\approx p_d(\mathbf{x})=\sum _{i=1}^r{\alpha }_i{\Psi }_i(\mathbf{x}),$$

$$\int {\psi }_k^{(j)}(x_j){\psi }_l^{(j)}(x_j){\rho }_j(x_j)d{x}_j={\delta }_{kl},$$

$${\Psi }_i(\mathbf{x})=\prod _{j=1}^d{\psi }_{i_k}^{(j)}(x_j),$$

$$\mathbf{i}=(i_i,...,i_d)\in {\mathbb{N}}_0^d.$$

$$p_d(\mathbf{x})=\sum _{i=1}^r{\alpha }_{\mathbf{i}}{\Psi }_{\mathbf{i}}(\mathbf{x}),{\Psi }_{\mathbf{i}}(\mathbf{x})=\prod _{j=1}^d{\psi }_{i_k}^{(j)}(x_j).$$

$$\left\{\mathbf{i}\in {\mathbb{N}}_0^d|\sum _{k=1}^d{i}_k\le p\right\},$$

$$\left\{\mathbf{i}\in {\mathbb{N}}_0^d|i_k\le p\text{for all}1\le k\le d\right\}.$$

$$\text{nz}((0,1,2,0,4))=\{2,3,5\}.$$

$$\sigma =\mathbb{E}\left[(p_d(\mathbf{x})-\mathbb{E}\left[p_d(\mathbf{x})\right]{)}^2\right]={\int }_{\mathcal{D}}{\left(\sum _{i=2}^r{\alpha }_{\mathbf{i}}{\Psi }_{\mathbf{i}}(\mathbf{x})\right)}^2\rho (\mathbf{x})d\mathbf{x}.$$

$$f(\mathbf{x})=\sum _{S\subseteq [d]}{f}_S({\mathbf{x}}_S)$$

$$f_S({\mathbf{x}}_S)=\mathbb{E}[f(\mathbf{x})|{\mathbf{x}}_S]-\sum _{T\subset S}{f}_T({\mathbf{x}}_T).$$

$${\sigma }_S=\frac{{\sum }_{\text{nz}(\mathbf{i})=S}{\alpha }_{\mathbf{i}}^2}{{\sum }_{\mathbf{i}}{\alpha }_{\mathbf{i}}^2}.$$

$$\widehat{{\sigma }_i}=\sum _{i\in S}{\sigma }_S.$$

$$\sigma {\sigma }_S=\left({\int }_{\mathcal{D}}{\int }_{\mathcal{D}}f(\mathbf{x})f({\mathbf{x}}_S;{\mathbf{x}}_{-S}^{(1)})d\mathbf{x}d{\mathbf{x}}^{(1)}-{\mu }^2\right),$$

$$t_S^{\prime }={\int }_{\mathcal{D}}{\int }_{\mathcal{D}}{\int }_{\mathcal{D}}f(\mathbf{x})f({\mathbf{x}}_S;{\mathbf{x}}_{-S}^{(1)})f({\mathbf{x}}_S;{\mathbf{x}}_{-S}^{(2)})d\mathbf{x}d{\mathbf{x}}^{(1)}d{\mathbf{x}}^{(2)}-{\mu }^3.$$

$$t_S^{\prime }+{\mu }^3=\mathbb{E}\left[\mathbb{E}[f(\mathbf{x})|{\mathbf{x}}_S{]}^3\right],$$

$$\sigma {\sigma }_S+{\mu }^2=\mathbb{E}\left[\mathbb{E}[f(\mathbf{x})|{\mathbf{x}}_S{]}^2\right].$$

$$S^1=\{p:\text{nz}(\mathbf{p})=S\},$$

$$S^2=\{(p,q):\text{nz}(\mathbf{p})\cup \text{nz}(\mathbf{q})=S\},$$

Full text

In this paper, we build on these well-known ideas by examining two different aspects of this paradigm. First, we study whether sensitivity indices can be computed efficiently if one leverages a polynomial ridge approximation—a polynomial fit over a subspace. Given the assumption of anisotropy in the dependence of a function, we show that sensitivity indices can be computed with a reduced number of model evaluations. Second, we discuss methods for evaluating sensitivities when constrained near output extrema. Methods based on the analysis of skewness are reviewed and a novel type of indices based on Monte Carlo filtering (MCF)—extremum Sobol’ indices—is proposed. We combine these two ideas by showing that these indices can be computed efficiently with ridge approximations, and explore the relationship between MCF-based indices and skewness-based indices empirically.

1 Introduction and motivation

Across all data-centric disciplines, sensitivity analysis is an important task. The phrase “which of my parameters is the most important?” can be heard across shop floors, office spaces and board rooms. In global sensitivity analysis, the measure of importance is typically the conditional variance, leading to the familiar Sobol’ indices [39, 40], the related total indices [16] and the Fourier amplitude sensitivity test (FAST) indices [6]. However, these are not the only metrics available: over the years, other metrics for sensitivity analysis have been proposed, including those based on conditional skewness [45] and kurtosis [7]. Beyond this, Hooker [17] and Liu and Owen[24] introduce the notion of a superset importance, which is the sum of all Sobol’ indices based on a superset of an input set. Techniques based on partial derivative evaluations of a function have also been extensively studied [3, 25], and the relationship between these elementary effects and total Sobol’ indices have been investigated. For example, in [21], a theoretical link between elementary effects and total Sobol’ indices is established. More recently, Razavi et al. [34, 33] provide a unifying framework for the analysis of sensitivities at different scales within the input domain via a variogram-based approach. In this approach, elementary effects and Sobol’ indices are recovered as special cases in the limit of small and large scales respectively. By analyzing sensitivities across the entire scale in the input space, structural features of the response surface, such as smoothness and multi-modality, are revealed.

In this paper, we are interested in the analysis of sensitivity near output extrema: which input variables—under what distributions—predominantly drive the output response near its highest/lowest levels? Output extrema are regions within the output space that are of special importance in reliability engineering; systems are often designed to operate under near-optimal conditions to maximize performance metrics. Extremum sensitivity analysis is useful for informing the construction of tailored surrogate models that aid optimization, by isolating critical variables that most influence the output near extrema. Moreover, the quantification of uncertainties at output extrema can help engineers understand the occurrence of extreme events, and how to mitigate the associated risks. The goal of this paper is to develop indices more localized than global sensitivity metrics, automatically adapted to regions of interest within the input domain corresponding to output extrema.

In the literature, papers directly addressing extremum sensitivity analysis are scarce. Recently, there has been work on higher-order indices utilizing skewness and kurtosis, directed at the analysis of the tails of the output distribution. In Owen et al. [29], the authors mention that a positively skewed distribution will attain more extreme positive values than a symmetric distribution with the same variance; distributions with a large kurtosis will also attain greater extremes on both sides [28]. Geraci et al. [11], who introduce polynomial chaos-based ideas for computing skewness- and kurtosis-based indices, find applications of these indices in more accurate approximation of the probability density function (PDF) of the output response, especially near output extrema. However, neither work directly addresses the problem of sensitivity analysis near output extrema. In this work, a novel type of sensitivity indices called extremum Sobol’ indices is proposed. As the name suggests, the aim is to compute Sobol’ indices of functions when the input distribution is restricted to regions in the domain leading to output extremum. This conditional distribution is modelled via a rejection sampling heuristic also known as Monte Carlo filtering (MCF).

Central to our investigation is the tractable computation of various sensitivity metrics. It is well-known that computing sensitivity indices—particularly in high-dimensional problems—can be prohibitive owing to the excessively large number of evaluations required. While access to a function’s gradients can abate this curse of dimensionality, our focus will be on the general case where access to a function’s gradients (or adjoint) is not available. Our approach is to use orthogonal polynomial least squares approximations of the uncertain function, commonly known as generalized polynomial chaos expansions (PCE) [12, 48, 44]. The application of PCE for global sensitivity analysis has been relatively well-studied and prior work can be found in [41, 1, 42, 37] and the references therein. The use of polynomial approximations implies the assumption that we are restricting our study to functions that can be well-approximated by globally smooth polynomials.

The remainder of this paper is structured as follows. In Section 2 we interpret the global polynomial approximation problem as that of computing its coefficients. For polynomials of a few dimensions, this can be done without assuming anisotropy; for higher-dimensional approximations, on the other hand, anisotropy can be used to circumvent the curse of dimensionality. More specifically, we construct polynomial approximations over subspaces—composed by a few linear combinations of the input variables. From the coefficients in the projected space, we design a polynomial surrogate for determining the coefficients in the full space, which lead to the familiar Sobol’ indices. In Section 3, we describe and provide algorithms for computing extremum Sobol’ indices using polynomials and their ridges. The methods and utility of extremum sensitivities are demonstrated through numerical examples in Section 4.

2 Sensitivity analysis with polynomial ridge approximations

Prior to detailing techniques for computing coefficients of the aforementioned polynomials, it will be worthwhile to make more precise our notion of anisotropy. Let us consider a polynomial

[TABLE]

For the case where $n=d$ and where $\bm{M}=\bm{I}$ (the identity matrix), we have a polynomial approximation in the full $d$-dimensional space. When $n<d$ and when the $n$ columns of $\bm{M}$ are linearly independent, they span an $n$-dimensional subspace in $\mathbb{R}^{d}$. In this case, the function $p$ is referred to as a generalized ridge function [32]. Depending on what the weights in $\bm{M}$ are, the polynomial will be anisotropic with the larger weights prescribing linear combinations along which the polynomial varies the most. For notational clarity, we add a subscript below $p$ to denote the dimension over which the polynomial is constructed: $p_{d}$ denotes a polynomial in the full space, while $p_{n}$ denotes a polynomial over a subspace. From here on, we will assume that $n\ll d$.

2.1 Computing polynomial coefficients using ridge approximations

Now, let us assume the existence of a smooth function $f:\mathbb{R}^{d}\rightarrow\mathbb{R}$ with input variables $\mathbf{x}$. Our goal is to compute sensitivities with respect to individual or groups of input variables. Provided one has access to $N$ input/output realizations

[TABLE]

and one can express $f$ as a polynomial (albeit approximately), we can state

[TABLE]

where $\bm{A}_{d}\in\mathbb{R}^{N\times r}$ is a polynomial design matrix in $d$ dimensions. Writing the expansion of $p_{d}$ in terms of its $r$ basis terms,

[TABLE]

we can identify $\bm{A}_{d}$ to be

[TABLE]

In other words, $\bm{A}_{d}$ is a Vandermonde-type matrix. The number of rows $N$ corresponds to the number of model evaluations, while the number of columns $r$ corresponds to the number of unknown polynomial coefficients. Consider two scenarios for the size of $N$ vs $r$:

$N\geq r$: This leads to a polynomial least squares formulation where one tries to solve

[TABLE]

The exponential scaling of $r$ with respect to $d$ requires $N$ to be chosen to match this scaling, implying that this approach is not computationally practical for moderate- to high-dimensional problems. In light of this, new strategies for selecting the placement of sample points have been an ongoing research topic, e.g. see [36, 37]. Recently, there has also been research into leveraging multi-fidelity to reduce the number of required evaluation points [30, 18].

• 2.

$N<r$: To solve a linear system with more unknowns than equations, one has to impose some form of regularization. Assuming that the solution is sparse on a polynomial basis, Blatman and Sudret [2] use least angle regression [10] to identify a subset of dominant terms in the basis. The assumption of compressibility leads to formulations based on compressed sensing [9], where strategies including subspace pursuit [8] and basis pursuit denoising [43, 13] have been proposed. The latter is an $\ell_{1}$-minimization strategy that uses second order cone programming to solve

[TABLE]

where $\epsilon$ is a small positive constant.

So, why the need to use a global polynomial approximation? As alluded to already, accurate polynomial approximations of $f$ permit us to leverage prior work [41] on computing global variance-based sensitivities using only the coefficient values. In this paper, we will also demonstrate how polynomial approximations are useful in computing Sobol’ indices tailored to output extrema.

But what if $f$, in addition to being smooth, was also anisotropic? The idea of developing global sensitivity metrics for ridge approximations has been studied previously by Constantine and Diaz [5]. Using ideas from active subspaces [4], the authors define activity scores that are based on the eigendecomposition of the average outer product of the gradient $\nabla f$ with itself in the input space. They draw comparisons between these activity scores and total Sobol’ indices (see Theorem 4.1 in [5]). While their work defines a new type of sensitivity index tailored to ridge approximations, our work builds upon the idea of ridge approximations with polynomials to propose a method for calculating moment-based sensitivity indices, such as the Sobol’ indices.

Consider the polynomial ridge approximation

[TABLE]

where $\bm{A}_{n}\in\mathbb{R}^{N\times q}$ is another polynomial design matrix, but now in $\mathbb{R}^{n}$. Since typically $q\ll r$, we will assume that $N\geq q$, implying that we have at least as many model evaluations as unknown coefficients in $\bm{\beta}$. Thus, computing $\bm{\beta}$ using regular least squares should not be a problem111We are assuming that the $N$ points are selected such that $\bm{A}_{n}$ has a relatively low condition number.. Note that the choice of the polynomial basis in the projected space is usually not critical. Provided the dimensionality of the projected space is small, the amount of training data required for finding the subspace is often more than sufficient for regression in the projected space. If this assumption is not true, a polynomial basis matching the distribution in the projected space can be constructed (see Jakeman et al. [19, sec. 5.6.2]) to mitigate losses of accuracy at this stage. In addition, in our case, sensitivity indices and moments are not quantified in the projected space, so the basis does not need to be matched to the distribution induced in the projected space (see Section 2.2).

The trouble, however, is that $\bm{\beta}$ are the coefficients of $p_{n}$, and we need access to the coefficients $\bm{\alpha}$ associated with $p_{d}$ in order to estimate global sensitivities. The key is to notice that $p_{n}(\bm{M}^{T}\mathbf{x})$ is simply a rearrangement of coefficients from $p_{d}(\mathbf{x})$. In other words, we can convert the coefficients of $p_{n}$ to those of $p_{d}$ in a lossless manner. Sampling $N\geq r$ points to evaluate the polynomial design matrices $\bm{A}_{d}$ and $\bm{A}_{n}$, the linear system

[TABLE]

can be solved to recover the full space coefficients $\bm{\alpha}$. Using this method, the accuracy of the estimated polynomial coefficients is now dependent on the accuracy of estimating the dimension-reducing subspace spanned by $\bm{M}$ and the error of a ridge approximation due to variations of $f$ outside the range of $\bm{M}$. Assuming the presence of strong anisotropy, the subspace can be computed with significantly fewer data than a full polynomial approximation using gradient-free techniques including polynomial variable projection [15] and minimum average variance estimation (MAVE) [47]. In the case that the input dimension is so large that $r$ is prohibitively large for solving this linear system, subset selection schemes for truncating the full space polynomial basis—such as that proposed by Blatman and Sudret [2]—can be used when determining the full space coefficients. This reduces the number of coefficients to be determined and prunes the basis of unnecessary terms, but the rearrangement step is no longer strictly lossless.

2.2 From polynomial coefficients to Sobol’ indices

Now, assume that there is uncertainty associated with the input variables. That is, we consider the inputs to be a set of random variables, whose distribution can be described by the PDF $\rho(\mathbf{x})$ with $\mathbf{x}\in\mathcal{D}\subset\mathbb{R}^{d}$ where $\mathcal{D}$ is the support of the input PDF. In this section, it is assumed that the set of inputs is independently distributed, with marginal probability densities $\rho_{j}(x_{j})$ for $j=1,...,d$ such that

[TABLE]

The Sobol’ indices are then proportions of the output variance $\sigma=\text{Var}[f(\mathbf{x})]$ attributed to individual or groups of input variables [39], defined for a set of variables $S$ as

[TABLE]

where $\mathbf{x}_{S}$ contains the variables in $S$. Directly evaluating this expression for a function via Monte Carlo requires a large number of function evaluations that scales poorly with input dimension. Below, we follow the work by Sudret [41] and show that Sobol’ indices can be computed via a PCE

[TABLE]

given a suitable orthogonal polynomial basis $(\Psi_{i})$.

Let $(\psi^{(j)}_{k})_{k\in\mathbb{N}}$ be an orthogonal basis of univariate polynomials with respect to $\rho_{j}(x_{j})$, where $k$ indicates the degree of the polynomial. The orthogonality of the basis implies

[TABLE]

where $\delta_{kl}$ is the Kronecker delta which equals to 1 if $k=l$ and 0 otherwise. Now, each multivariate basis function in 12 $\Psi_{i}$ is a product of univariate orthogonal basis polynomials in each dimension, i.e.

[TABLE]

for all $i$ where $i_{k}$ is the degree of the polynomial in the $k$-th dimension. Thus, define a multi-index that associates each multivariate basis polynomial with the degrees of its constituent univariate basis polynomials,

[TABLE]

The polynomial series 12 can now be rewritten using multi-indices,

[TABLE]

The multi-index enables us to conveniently specify the polynomials included in our basis. The set of included multi-indices is called the index set. For instance, an isotropic index set of maximum total order $p$ is specified by

[TABLE]

and a tensor grid of order $p$ is specified by

[TABLE]

Furthermore, define a function $\text{nz}(\cdot)$ that returns the positions of the non-zero indices in a multi-index. For instance,

[TABLE]

Leveraging the notation introduced so far, we can express the decomposition of the (approximate) total output variance $\sigma=\text{Var}[p_{d}(\mathbf{x})]$ via the following equation

[TABLE]

Owing to orthonormality, the polynomial series expression 16 allows us to compute the Sobol’ indices straightforwardly given the polynomial coefficients $\bm{\alpha}$ by comparison with the Sobol’ decomposition

[TABLE]

where $[d]=\{1,2,...,d\}$, $\mathbf{x}_{S}$ selects all coordinates of $\mathbf{x}$ belonging to the subset of indices $S$, and $f_{S}(\mathbf{x}_{S})$ is recursively defined as

[TABLE]

Using this decomposition, Sudret shows that the Sobol’ index for a set of variables indexed by $S$ is given by

[TABLE]

A measure of the importance of an input variable is the total Sobol’ index [40], which is the sum of all Sobol’ indices that contain the variable concerned. This index factors in all orders of interaction within the function that involves the concerned variable. For instance, for variable $x_{i}$, the corresponding total Sobol’ index is

[TABLE]

3 Extremum sensitivity analysis

In the previous section, we described a method to use coefficients obtained from a polynomial ridge approximation to calculate Sobol’ indices. In the following section, the second key contribution of this paper—the proposal of polynomial-based sensitivity indices focusing on inference near output extrema—will be addressed.

3.1 Skewness-based indices

First, we review two works in literature related to global sensitivity analysis based on the decomposition of skewness. Assuming that the input domain is a unit hypercube $\mathcal{D}=[0,1]^{d}$ endowed with uniform marginal distributions, Owen et al. [29] define higher-order Sobol’ indices by generalizing the Sobol’ identity, which is based on the functional ANOVA (Sobol’) decomposition (21). This identity states that [40, Thm. 2]

[TABLE]

where the notation $\mathbf{x}_{S};\mathbf{x}^{(1)}_{-S}$ denotes a point $\mathbf{y}$ in $\mathcal{D}$ where $y_{j}=x_{j}$ for $j\in S$ and $y_{j}=x^{(1)}_{j}$ for $j\notin S$, and $\mu$ denotes the mean of the output. The generalization for calculating the (unnormalized) skewness index for a subset of variables $S$ reads

[TABLE]

Also in [29], it is shown that

[TABLE]

which parallels the definition of Sobol’ indices where

[TABLE]

The authors in [29] remark that skewness indices may provide indication that certain variables are important for attaining extreme values of the output, with the sign of the skewness indices indicating whether the variables contribute significantly to output maximum (positive skewness index) or minimum (negative skewness index). However, they do not provide further theoretical justification for this claim.

The formulation described by Owen et al. lends straightforwardly to computation by Monte Carlo samples, but the high dimensionality of the integrals implies that many samples are required. In Geraci et al. [11], skewness indices are defined by evaluating the third central moment of the Sobol’ decomposition (21) and expanding the sum, resulting in indices that are different from those defined in [29]. For brevity, we omit the full exposition and refer the reader to their work for further details. The authors provide an approach for evaluating these indices via polynomial chaos, by expressing the indices in terms of polynomial coefficients. However, unlike the case for variance-based indices, the resultant indices cannot be expressed in a form as simple as 23. Here, we cite the full result from [11] for the skewness index of a set of variables indexed by $S$. Letting $\mathbf{p},\mathbf{q},\mathbf{r}$ be the corresponding multi-indices to the scalar indices $p,q,r$ respectively, we define the following sets

[TABLE]

[TABLE]

[TABLE]

It can be shown that the skewness index for $S$ is

[TABLE]

The total skewness index for $x_{i}$ is defined analogously to 24

[TABLE]

In the same work, the authors also establish the importance of these indices in modelling the tails of the output PDF, but make no explicit reference to sensitivity analysis near output extrema.

3.2 Extremum Sobol’ indices

Although skewness indices pertain to the analysis of the tails of the output distribution, limitations exist with respect to the use of these indices for extremum sensitivity analysis. For a simple example, consider the following function

[TABLE]

By inspection, it is clear that $x_{2}$ is responsible for none of the (small but non-zero) output skewness, but it is more important than $x_{1}$ throughout the entire input domain, including regions of output extrema.

In this paper, we propose a novel set of variance-based metrics that quantify the relative input sensitivities near regions in the domain yielding extreme output values, called extremum Sobol’ indices. The aim is to analyze the variance decomposition of the function conditioned near its extrema, producing indices that have a clear interpretation for extremum sensitivity analysis. To characterize the input distribution corresponding to output conditioning, a heuristic related to Monte Carlo filtering (MCF) [35, p. 39] is used. In MCF, the input space is partitioned into two regions—one where the corresponding outputs fall in a range of interest $B$, and one where this does not happen, $\bar{B}$. The importance of a variable $x_{i}$ is determined by the difference between the distributions of $x_{i}$ conditioned within $B$ and its complement $\bar{B}$ respectively, where the distributions are inferred using Monte Carlo. In our method, MCF is used to isolate regions containing points leading to output maximum and minimum, and the Sobol’ indices of the quantity of interest are evaluated under the conditional distributions restricted to these regions. Sobol’ indices resulting from points leading to output maximum are called top Sobol’ indices, and those from points leading to output minimum are called bottom Sobol’ indices.

Computation of extremum Sobol’ indices involves the following steps.

Global polynomial approximation: Fit a global polynomial approximation to the model by input/output pairs $\{\mathbf{x},f(\mathbf{x})\}$ via sampling from $\mathcal{D}$ according to $\rho(\mathbf{x})$. Note that this global approximation can take the form of a ridge approximation. 2. 2.

Rejection sampling: Generate another $N$ random input points in $\mathcal{D}$ according to $\rho(\mathbf{x})$. Estimate the maxima and minima by pushing these samples through the polynomial approximation. Isolate input-polynomial output $\left\{\mathbf{x},p(\mathbf{x})\right\}$ pairs that fall within the top 5% and bottom 5% of the maxima and minima respectively. 3. 3.

Marginals and correlation determination: Approximate the marginal distributions for each parameter for the two groups using input/output pairs. Fit a tailored bandwidth kernel density estimate to these marginals. Additionally, evaluate the empirical correlation matrix for these two groups. 4. 4.

Extremum Sobol’ indices: Construct polynomial chaos approximations over the groups, using additional input/output pairs $\{\mathbf{x},f(\mathbf{x})\}$ sampled from the correlated distribution over the respective groups. Then, use Algorithm 1 to compute the extremum Sobol’ indices for each group.

For the final step, it is important to note that the extremum points resulting from MCF are in general distributed with a correlated measure, even if the original input measure is independent. In this case, the output variance cannot be decomposed in a similar manner to the Sobol’ decomposition, and the indices computed from the definition (11) mixes the contribution of variables outside of the variables in $S$; the sum of all Sobol’ indices thus does not result in unity. Therefore, we compute sensitivities using the formulation proposed by Li et al. [22] which involves a covariance decomposition, for which Navarro et al. [27] expanded upon with an algorithm based on polynomial chaos expansions. Algorithm 1 in A summarizes the steps of this computation given extreme points isolated by MCF.

To calculate the extremum Sobol’ indices using a polynomial chaos approximation, a polynomial basis orthogonal to the correlated measure determined in step 4 above is required. There are different methods to construct orthogonal polynomial bases in correlated spaces. For instance, one can map the dependent variables to an independent space via the Nataf transform [23]. In Algorithm 1, we use the Gram-Schmidt method [46, 19] to generate evaluations on a polynomial basis orthogonal to the correlated measure. We start with the polynomial design matrix $\bm{A}_{u}$ as in equation (5), with the $N$ points of evaluation $\mathbf{x}_{i}$ sampled from the correlated measure. The basis polynomials here can be taken as ones orthogonal to the full input space measure $\rho(\mathbf{x})$. The sampling step can be performed via the Nataf transform. Then, the QR decomposition

[TABLE]

is calculated, where $\bm{Q}_{u}$ has orthogonal columns, and $\bm{R}$ is upper triangular. Because of the distribution of the sample points, this step is a Monte Carlo approximation to the orthogonalization of the polynomial basis with respect to the inner product

[TABLE]

where $\omega$ is the correlated (extremum) measure and $\Psi_{1}$ and $\Psi_{2}$ are two polynomial basis functions. Polynomial coefficients can be solved on this basis using least squares and Sobol’ indices quantified via the coefficients. However, note that the Gram-Schmidt process produces basis polynomials that mix contributions from different variables. Therefore, the Sobol’ indices need to be calculated via Monte Carlo on the basis polynomials instead of being simply read out from (23). For this algorithm, we remark that additional function evaluations are required only for solving the coefficients in the new polynomial basis; all Monte Carlo steps are performed on evaluations of the polynomial basis terms.

4 Numerical examples

In the following section, polynomial ridge approximations and extremum sensitivity analysis are applied on several numerical case studies. We aim to highlight the following points:

Demonstrate the estimation of Sobol’ indices using polynomial ridge functions via the least squares formulation described in Section 2. 2. 2.

Discuss empirical parallels between skewness indices (Section 3.1) and extremum Sobol’ indices (Section 3.2) at capturing important features at top and bottom output values. 3. 3.

Compare the computational cost of our framework with Monte Carlo and other polynomial chaos methods such as adaptive LARS where applicable.

4.1 Ridge approximation of an analytical function

In this example, we study the estimation of Sobol’ indices of the function $f:[-1,1]^{6}\rightarrow\mathbb{R}$, where

[TABLE]

with

[TABLE]

[TABLE]

and the orthogonal columns of $\bm{U},\bm{W}$ span, respectively, the column and left null space of

[TABLE]

Here, $s$ is a parameter that modulates the deviation of $f(\mathbf{x})$ from an exact ridge function within the column space of $\bm{U}$—i.e. it allows the function to vary slightly outside of the plane spanned by the columns of $\bm{U}$. We evaluate the sensitivity indices from noisy samples $\{\mathbf{x}_{i},f(\mathbf{x}_{i})+\epsilon_{i}\}_{i=1}^{N}$, where the sample inputs $\mathbf{x}_{i}$ are drawn from the uniform distribution over the input domain, $\rho(\mathbf{x})=1/2^{6}$, and $\epsilon_{i}$ represents additive Gaussian random noise with a standard deviation of $10^{-4}$. Five methods are compared for this task:

Polynomial ridge approximation: Using MAVE [47]—a gradient-free dimension reduction method—we estimate a dimension-reducing subspace for the function of interest. Then, we fit a low-dimensional quadratic polynomial ridge over the reduced space $p_{n}$ with $n=2$, whose coefficients we use to calculate the coefficients in the full space via 9. Since $n$ is small, we can afford to use a tensor grid. A Legendre polynomial basis is used in the projected space; since we do not need to quantify moments, the choice of basis is not critical here. Note that in this example, the choice of $n=2$ is prescribed according to the known function, but in practice for an unknown model, methods such as cross validation (or visualization with projection plots) may be used for determining the suitable dimensionality.

• 2.

Quasi-Monte Carlo: We compute the integrals defining the Sobol’ indices via a Sobol’ sequence, as specified in [40]. We adapt the open-source Python library SALib [14] for this method.

• 3.

Adaptive least angle regression (LARS): Based on an approach by Blatman and Sudret [2] (see C), we evaluate the Sobol’ indices via orthogonal polynomial regression over a pruned basis, which is a subset of an isotropic total order basis of maximum degree 4 over the full space. The open-source Python library scikit-learn [31] is used for LARS.

• 4.

Fitting a polynomial on the full input space via least squares with a total order basis of maximum order 3.

• 5.

Fitting a polynomial on the full input space via least squares with a total order basis of maximum order 4.

Note that the final two approaches place a lower limit on the number of samples required (84 and 210 respectively) as explained in section 2.1. The two largest Sobol’ indices $\sigma_{13}$ and $\sigma_{24}$ are computed using the methods above. To benchmark these results, we also compute the quantities using integration via Gauss quadrature on a degree 4 tensor grid of Legendre polynomials with $(4+1)^{6}=15625$ points, the results of which we treat as the truth. The fact that $f(\mathbf{x})$ is a polynomial allows this integration to be exact. Figure 1 shows the results of this comparison for $s=0$ and $s=0.1$, respectively.

From the results, it is clear that polynomial-based techniques offer a significant advantage over quasi-Monte Carlo integration. This shows that sensitivity analysis of a function that is smooth and well approximated by a polynomial can be greatly facilitated by our approach. Additionally, the fact that this function is (approximately) a ridge function gives the ridge approximation approach an advantage over other sparse methods, including adaptive least angle regression. In general, ridge functions need not have sparse coefficients unless the ridge structure aligns well with the canonical directions in the input space. The case of $s=0.1$ adds variation outside of the two-dimensional subspace spanned by the columns of $\bm{U}$, which reduces the accuracy of the ridge approximation. Nonetheless, the ridge approximation approach fares well against other sparse approaches. However, if we are allowed to evaluate the function more than 210 times, regression on the full space using a total order grid is superior.

4.2 Extremum sensitivity analysis of an analytical function

Consider the additive function $f:[0,1]^{2}\rightarrow\mathbb{R}$ where

[TABLE]

We calculate total sensitivity indices for this function via polynomial methods, which are exact in this case because the function is a polynomial222However, Monte Carlo is still used to calculate the extremum Sobol’ indices after polynomials are fitted near output extrema. Sufficient samples are used to ensure that the uncertainty is small.. Sobol’ indices, skewness indices and extremum Sobol’ indices are plotted for this function in Figure 2. In all bar chart figures, the bar heights show the mean and red error bars show $\pm 1$ standard deviation over 30 trials. It can be seen from the middle plot that both variables have similar Sobol’ indices, but there is a clear difference for the bottom and top Sobol’ indices. For low output values, $f_{2}$ is nearly flat, implying that it is not influential to the function value; at high outputs, $f_{1}$ is relatively flat, implying that $f_{2}$ should be the more important variable. These pieces of information are revealed from the extremum Sobol’ indices.

From the skewness indices, it can be seen that $x_{1}$ is responsible for skewing the output distribution to the right, and $x_{2}$ to the left. In this example, it can be observed that $x_{1}$, with a negative skewness index, is important near the output minimum; for $x_{2}$, it is important near the output maximum, and has a positive skewness index. This agrees with the interpretation of skewness indices by Owen et al. [29] mentioned in Section 3.1, although it should be noted that their definition of skewness indices is different from the ones computed here. However, whether the sign of skewness indices are indicators of extremum sensitivities in general, and its precise connection with extremum Sobol’ indices remain open topics.

4.3 Borehole function

In this case study, we apply methods of extremum sensitivity analysis on the borehole function, which models the steady-state water flow rate for a hypothetical borehole geometry [26]. Our objective here is to identify which variable is responsible for the greatest variation over extreme values of the steady-state water flow. The model is specified by eight input parameters assumed to be independent. Their distributions are given in Table 1. The water flow rate depends on the parameters as

[TABLE]

Treating the water flow rate as the quantity of interest, we compute the total Sobol’ and skewness indices with respect to the input variables. Three methods of computation are compared:

(Poly) Using an orthogonal polynomial least squares approximation on an isotropic basis with maximum total order of 3 (with a cardinality of 165), trained from 300 input/output pairs. 2. 2.

(Ridge) Using a polynomial ridge approximation of degree 3 in a 4-dimensional subspace333In this example and the next, the subspace dimension is found by trial-and-error, and the lowest value that results in acceptable performance is selected. In practice, methods such as cross validation (or visualization with projection plots) can need to be used for determining the suitable dimensionality, but we omit a detailed study here for brevity., whose coefficients are subsequently converted for a polynomial on an isotropic basis with maximum total order of 3. The subspace is found via polynomial variable projection [15] using 150 input/output pairs. In this and the next numerical examples, Legendre bases on tensor grids are used for fitting polynomials in projected spaces. 3. 3.

(MC) Using $10^{5}$ model evaluations to calculate the indices via the Monte Carlo method by Sobol’ [40, Thm. 2].

The training set $\{\mathbf{x}_{i},f(\mathbf{x}_{i})\}_{i=1}^{N}$ is formed by sampling the inputs independently according to the corresponding marginal input distributions, and evaluating the corresponding output values. In Figure 3, results for the computations are shown. The third method is omitted for the skewness indices because no equivalent Monte Carlo formulation is proposed for skewness indices. The algorithm for computing skewness indices via polynomial chaos is explained in B. From this figure, it can be seen that the Sobol’ indices computed using a polynomial approximation matches well with the Monte Carlo estimates, and using a ridge approximation in lieu of a full space polynomial reduces the amount of training data required while having minor effects on the accuracy of the computed indices. The choice of polynomial degrees for this and the next example is confirmed by verification results found in D.

Next, we compute the extremum Sobol’ indices using the procedure described in Section 3.2. The total top and bottom Sobol’ indices are calculated using MCF to select 5% of the input points yielding maximum and minimum outputs respectively out of $10^{5}$ points. Three variants of the computation algorithm are tested differing in the form of the model used to generate the pool of samples for MCF at step 2—namely

Using $10^{5}$ function evaluations directly (instead of a polynomial approximation as specified in step 1); 2. 2.

Using a global polynomial model at total degree 4 (cardinality 495) fitted with 700 observations; and 3. 3.

Using a ridge polynomial at total degree 4 in a 4-dimensional subspace fitted with 400 observations.

Note that the polynomial degree is increased by one compared to the approximation models used to calculate global sensitivity indices previously in this section. This is to ensure better accuracy at the tails of the output. For all methods, an additional 1400 function evaluations are used to fit a degree 5 polynomial at the extrema.

The results over 30 trials are shown in Figure 4. It can be seen that near high outputs values, $r_{w}$ has significantly increased importance compared to that predicted by full space Sobol’ indices, along with somewhat increased importance for $H_{u}$, $H_{l}$ and $L$. For low outputs, the importance ranking is similar to that evaluated from full space Sobol’ indices. In both cases, it can be seen that no significant loss of accuracy is incurred by using either polynomial or ridge approximations in place of function evaluations for MCF. In Figures 5 and 6, the marginal and pairwise distributions obtained from MCF are plotted, showing some correlation between $r_{w}$ and $K_{w}$ which is not present in the original full space input distribution. From the skewness indices, we see the increased importance of $r_{w}$, $H_{u}$, $H_{l}$ and $L$ with positive values. This corroborates with the interpretation of skewness indices mentioned in Section 4.2. The agreement is only qualitative however, with the relative ranking of $r_{w}$ and $K_{w}$ not in agreement. Again, it is not clear whether the observations can hold generally.

4.4 Piston cycle time function

Consider a model of a piston adapted from [20], which describes the cycle time of the piston as a function of seven parameters. These parameters are shown in Table 2. The inputs are independently and uniformly distributed over their support. As in the previous example, we are interested in obtaining the relative importance of each input variable near output extrema. The cycle time in seconds is

[TABLE]

where

[TABLE]

and

[TABLE]

The total Sobol’ and skewness indices are evaluated using the same three methods as in the previous example—Poly, Ridge and MC. For Poly, a degree 4 polynomial on a total order isotropic basis (cardinality 330) is fitted with 500 samples; for Ridge, a degree 4 polynomial over a 4-dimensional subspace is fitted with 300 samples and its coefficients. The results are shown in Figure 7, again showing that using ridge polynomial approximations incurs insignificant losses in accuracy.

The extremum Sobol’ indices are evaluated using the same approaches as in the previous example. For Poly, we use a degree 5 polynomial (cardinality 792) trained from 1000 samples; for Ridge, a degree 5 polynomial over a 4-dimensional subspace trained with 700 samples is used. A degree 4 polynomial in the extrema is fitted with an additional 500 training points for evaluating the indices. The results are shown in Figure 8. Comparing with full space Sobol’ indices, the bottom Sobol’ indices show an increased importance of $V$ in driving the output at low values. For high outputs, $k$ becomes important. Examining the total skewness indices, the increased importance at maximum can be seen for $k$, which has a large positive total skewness index. The information for $V$ is not found in the total skewness indices, but examining the first order skewness indices (Figure 9), the importance of $V$ at minimum is highlighted by the large negative value. From the top and bottom Sobol’ indices, we find that this variable is important throughout the domain including the output extrema at both ends, but it is especially important at the minimum. With skewness indices, it is necessary to examine different orders to draw this conclusion. Again, we note that in the above, empirical observations about skewness indices are described, but no attempt has been made for a unifying theoretical statement—an interesting topic that we invite further discussions about. Figures 10 and 11 show the pairwise plots of the extrema PDFs, showing correlations between $S$, $V$ and $k$ near extrema, and justifying the use of polynomials taking this into account.

5 Conclusions

In this paper we explore two related threads in sensitivity analysis. Through a polynomial-based approximation framework, we show the possibility of leveraging the ridge structure of a quantity of interest to mitigate the curse of dimensionality and reduce the number of observations required for sensitivity analysis. In addition, the use of polynomials and their ridges is demonstrated in the evaluation of a new type of sensitivity metric that quantifies sensitivities near output extrema—the extremum Sobol’ indices, reducing the cost of otherwise expensive Monte Carlo filtering operations. Since some empirical connections can be observed between extremum Sobol’ indices and skewness indices, we hope that this work will encourage discussion on the interpretation and utility of skewness indices in extremum sensitivity analysis.

Acknowledgements

This work was supported by the Cambridge Trust; Jesus College, Cambridge; the Alan Turing Institute and Rolls-Royce plc. We would also like to thank Art Owen for helpful discussions, and anonymous reviewers for their comments, which helped improve the manuscript.

Appendix A Computing extremum Sobol’ indices

The algorithm for computing extremum Sobol’ indices after obtaining extremum samples from MCF is shown in Algorithm 1. Note that Sobol’ indices beyond the first order need to be computed recursively by evaluating the associated quantities at lower orders (in line 15).

Appendix B Computing skewness indices

As described with (29) in Section 3, the method for calculating skewness indices via polynomial expansion coefficients is more involved than calculating the Sobol’ indices, because the integrals of third-order cross terms cannot be simplified easily with orthogonality. Hence, the expectation terms within each summand need to be evaluated using numerical quadrature. That is, we can approximate the expectation of a function $f(\mathbf{x})$ using Gauss quadrature

[TABLE]

where $(\mathbf{s}_{q},\nu_{q})_{q=1}^{M}$ are the quadrature points and weights. In Algorithm 2 we describe the algorithm implemented in Effective Quadratures [38]. In this algorithm, we have defined the normalized multi-index corresponding to the set of variables indexed by $S$

[TABLE]

where $\bar{j_{i}}=1$ if $i\in S$ and 0 otherwise, and $\odot$ denotes elementwise multiplication. In general, an $O(r^{3})$ loop is required for calculating a skewness index, and since $r$ can scale exponentially with the input dimension, this presents a significant computational burden. Our algorithm incorporates several methods to reduce the computational time, namely

Using a selection function similar to that proposed in [11] to avoid computation of terms which we know are zero;

• 2.

Scanning the index set in an $O(r)$ loop before computation to filter out indices that involve more variables than specified. This cuts the cost of computation significantly for indices involving a small number of variables.

• 3.

Using a formulation amenable to matrix operations to avoid explicit loops.

The code for this algorithm can be found at https://www.effective-quadratures.org.

Appendix C Adaptive least angle regression

In this section, the method of adaptive least angle regression (LARS) implemented for comparison in the numerical examples (Section 4) is briefly described. This method is based on [2], and the interested reader is referred to this work for further details. Define the predictor vector $\phi_{i}$ to be the evaluation of the $i$-th basis function on the training data set, i.e. the $i$-th column of the polynomial design matrix (5). The LARS algorithm [10] proceeds as follows:

Standardize all predictors and start with the residual $\delta=\mathbf{y}-\bar{\mathbf{y}}$. Initialize all coefficients $\bm{\alpha}=\mathbf{0}$. 2. 2.

Find the predictor $\phi_{i}$ most correlated with $\delta$. 3. 3.

Increase $\alpha_{i}$ until some other $\phi_{j}$ has as much correlation with $\mathbf{r}$ as $\phi_{i}$ does. 4. 4.

Increase both $\alpha_{i}$ and $\alpha_{j}$ towards their joint least squares coefficient on $\{\phi_{i},~{}\phi_{j}\}$ until another predictor $\phi_{k}$ has as much correlation with the current residual. 5. 5.

Repeat until min($r,~{}N$) predictors have entered, where $r$ is the number of basis functions and $N$ is the number of available data points.

The adaptive LARS procedure applies LARS for selecting a sparse basis, as follows:

Choose an initial candidate basis with maximum polynomial degree of $p=1$. 2. 2.

Apply LARS to the candidate basis and retain the terms with non-zero coefficients. 3. 3.

Calculate an estimate of the error of the model on this pruned basis via the leave-one-out estimate on the training data (see [2]). 4. 4.

Repeat for increasing $p$ and find the model with the least error at order $p^{*}$. 5. 5.

Using the basis found for $p^{*}$, fit a polynomial using least squares regression.

Note that this procedure is a special case of the method described in [2] in two respects: 1) A fixed data set is used instead of a sequential experimental design for fair comparison with other methods in this paper; and 2) An isotropic total order basis is used in all instances, consistent with other methods. For the results in this paper, the version of LARS implemented in the open-source Python package scikit-learn [31] is used.

Appendix D Additional verification results for the borehole and piston functions

In this section, the choices of polynomial degrees used for sensitivity analysis of the piston and borehole models are verified. First, Figure 12 verifies the accuracy of the extreme polynomials (polynomials fitted to extreme points from MCF) by plotting polynomial evaluations against true model evaluations for points near both extrema for both functions. The points shown are testing points that are not used in the fitting process. This figure shows that the polynomial degree chosen gives sufficient accuracy at the desired output range.

Second, the sensitivity indices for the borehole and piston functions are computed again with perturbations to the degrees of the polynomial approximations. Note that for extremum Sobol’ indices, the polynomial degree refers to that of the polynomial fitted to the extremum samples (according to the extremum distribution), not the polynomial used to obtain MCF samples. The latter is already compared with actual function evaluations which ensures its accuracy. Figures 13, 14, 15 and 16 show that perturbing the polynomial degrees produces negligible changes to the values of the indices, thus confirming that the degrees used in the main text are appropriate.