Derivative based global sensitivity analysis and its entropic link

The variance-based total effect measure accounts for the total contribution of an input to the output variation, and is often a preferred approach due to its intuitive interpretation and quantitative nature.

Let us denote $\mathbf{x}=(x_{1},x_{2},\dots,x_{d})$ as independent random input variables, and $y$ being the output of our computational model represented by a function $g$, such that $y=g(\mathbf{x})$.

The variance-based GSA decompose the output variance $V(Y)$ into conditional terms [2]:

where

and so on for the higher interaction terms. $V_{i}$ measures the first order effect variance and $V_{ij}$ for a second order effect variance, where their contributions to the unconditional model output variance can be quantified as $V_{i}/V(Y)$ and $V_{ij}/V(Y)$ respectively. Analogous formulas can be written for higher-order terms, enabling the analyst to quantify the higher-order interactions.

The total order sensitivity index is then defined as:

where $x_{\sim i}$ is the set of all inputs except $x_{i}$. The total order sensitivity index measures the total contribution of the input $x_{i}$ to the output variance, including its first order effect and its interactions of any order with other inputs.

The most direct computation of $S_{T_{i}}$ is via Monte Carlo (MC) estimation because it does not impose any assumption on the functional form of the response function. Various estimators for $S_{T_{i}}$ have been developed where the main difference is in the sampling design. A recent paper on comparing performance of several existing total-order sensitivity estimators can be found in [16]. Regardless of the specific estimator chosen, the computational cost for $S_{T_{i}}$ is a least $N(d+1)$, where $N$ is base sample number. Large values of $N$, normally in the order of thousands or tenths of thousands, is needed for more accurate estimate of $S_{T_{i}}$, and the computational cost has been noted as one of the main drawbacks of the variance-based GSA in [2].

The elementary effect method, also known as Morris’ method [2, 4], on the other hand is a simple and effective screening method for GSA. The Morris’ method is based on the construction of $r$ trajectories in the sampling space and the total cost of computation is $r(d+1)$. $r$ is typically between 10 and 50 [17], as compared to $N$ which is in the order of thousands for the estimation of $S_{T_{i}}$.

When the function $g(\cdot)$ is differentiable, functionals based on $g_{x_{i}}=\partial g/\partial x_{i}$ have been proposed to examine the global sensitivity with respect to the parameter $x_{i}$. The modified Morris sensitivity index $\mu^{*}$ [17] is an approximation of the following functional:

for $i=1,2,\dots,d$ and $f_{X}(\mathbf{x})$ is the joint probability density function (jPDF) of $\mathbf{x}$. It has been empirically demonstrated in [17] that for some practical problems $\mu^{*}$ is similar to the variance-based total index $S_{T_{i}}$.

Based on that observation, the Derivative-based Global Sensitivity Measure (DGSM):

has been proposed to be used as a proxy for $S_{T_{i}}$ [7, 8]. In particular, the total sensitivity variance $V(Y|\mathbf{X}_{\sim i})$ is upper bounded by DGSM via the following inequality:

which is developed in [8] for random inputs with log-concave distributions such as exponential, beta and gamma distributions. $c_{i}$ is called the Cheeger constant. Although the Cheeger constant has no analytical expressions for general log-concave measures, it can be estimated numerically using the distribution function of the uncertain input as shown in [8].

Divide Eq 5 by the output variance $V(Y)$ from both sides, we then have the upper bound $S_{T_{i}}\leq 4c_{i}^{2}\nu_{i}/V(Y)$ which can be used as a proxy for $S_{T_{i}}$. For Gaussian inputs with variance $\sigma_{i}^{2}$, the Cheeger constant $c_{i}=\sigma_{i}/2$ and the inequality in Eq 5 becomes $S_{T_{i}}\leq\sigma_{i}^{2}\nu_{i}/V(Y)$ as given in [7].

It has been pointed out in [18] that a major limitation of variance-based sensitivity indices is that they implicitly assume that output variance is a sensible measure of the output uncertainty, which might not always be the case. For instance, if the output distribution is multi-modal or if it is highly skewed, using variance as a proxy of uncertainty may lead to contradictory results.

To illustrate this point, we look at the simple function $y=x_{1}/x_{2}$. In this case, the two inputs both follow the chi-squared $\chi^{2}$ distribution with $x_{1}\sim\chi^{2}(10)$ and $x_{1}\sim\chi^{2}(13.978)$, and are assumed to be independent. This results in a positively skewed distribution of $Y$ with a heavy tail. This example has been used in [19] to show that variance-based sensitivity indices fail to properly rank the inputs of a model whose output has a highly skewed distribution.

We reproduce the sensitivity results of $S_{T_{i}}$ and $KL_{T_{i}}$ in Table 1 from [19]. In comparison, we have also computed the DGSM results, using numerical integrations as the derivatives are analytically available for this simple function. It can be seen from Table 1 that DGSM also provides misleading results in this case.

In addition, we also compare the results with the entropy-based total sensitivity index [20]:

where the average is calculated over all possible values of $\mathbf{X}_{\sim i}$ and $H$ is the differential entropy, e.g. $H(Y)=-\int f(y)\ln f(y)dy$. $\eta_{T_{i}}$ measures the remaining entropy of $Y$ if the true values of $\mathbf{X}_{\sim i}$ can be determined, in analogy to the variance-based total effect index $S_{T_{i}}$. Note that $\eta_{T_{i}}\leq 1$ because $H(Y|\mathbf{X}_{\sim i})\leq{H(Y)}$ [21].

$\eta_{T_{i}}$ has been estimated numerically based on the histogram method given in Appendix A, where $10^{7}$ samples are used. It can be seen from Table 1 that the entropy-based $\eta_{T_{i}}$ is able to effectively identify the higher influence of $X_{1}$ for this example. Note that different from $KL_{T_{i}}$ which is conditional on the value of $x_{i}$, for example $x_{i}$ are set at their mean values in Table 1, $\eta_{T_{i}}$ is an un-conditional sensitivity measure as all possible values of the inputs are averaged out.

We note in passing that analogously to the variance based sensitivity indices, a first order entropy index can also be defined as $\eta_{i}={H(Y)-H(Y|X_{i})}/{H(Y)}=I(X_{i},Y)/{H(Y)}$ [22]. $I(X_{i},Y)$ is the mutual information which measures how much knowing $X_{i}$ reduces uncertainty of $Y$ or vice versa. The index $\eta_{i}$ can thus be regarded as a measure of the excepted reduction in the entropy of of $Y$ by fixing $X_{i}$.

From the motivating example, it became clear that variance-based total sensitivity index might provide misleading sensitivity results, due to the fact that only the second moment is considered. This is especially the case for highly skewed or multi-modal distributions . This limitation is overcome by the entropy-based measures which are applicable independent of the underlying shape of the distribution.

However, the entropy-based indices have limited application in practice, mainly due to the heavy computational burden where the knowledge of conditional probability distributions are required. Both histogram and kernel based density estimation methods have computational challenges for entropy-based sensitivity indices as pointed out in [18]. For example, it was noted in [3] that although both variance-based and entropy-based sensitivity analysis take long computational time, the convergence for conditional entropy estimation is even slower.

Let $X$ be a continuous random variable with the probability density function (PDF) given by $f_{X}(x)$. $Y=g(X)$ can be regarded as a transformed variable and the transformed PDF $f_{Y}(y)$ can then be found via the Jacobian matrix for the calculation of its differential entropy $H(Y)=-\int f(y)\ln f(y)dy$.

For a general vector transformation $\mathbf{Y}=\mathbf{g}(\mathbf{X})$, the corresponding differential entropy are related as [23]:

where $\mathbb{J}$ is the Jacobian matrix with $\mathbb{J}_{ij}=\partial g_{i}/\partial x_{j}$ and $f_{X}(\mathbf{x})$ is the joint probability density function (jPDF) of $\mathbf{x}$. The above inequality becomes an equality if the transform is a bijection, i.e. an invertible transformation. Note that we have not assumed independent inputs for the inequality above.

As show in [23], Eq 8 can be proved substituting the transformed PDF $f(\mathbf{y})=f(\mathbf{x})/\det{\mathbb{J}}$ into the expression of $H(Y)$ and note there will be a reduction of entropy if the transformation is not one-to-one. Following this line of thought, it was noted in [24] that Eq 8 can also be seen as one version of the data processing inequality, where the transformation does not increase information.

Now set $y_{1}=g_{1}(\mathbf{x})=g(\mathbf{x})$ where $\mathbf{x}=(x_{1},x_{2},\dots,x_{d})$, and introduce dummy variables $y_{i}=g_{i}(\mathbf{x})=x_{i}$ with $i=2,\dots,d$. In this setting, the Jacobian matrix from the 2nd row onwards, i.e. $i\geq 2$, $\partial g_{i}/\partial x_{j}=1$ when $i=j$ and $\partial g_{i}/\partial x_{j}=0$ when $i\neq j$. Therefore, the Jacobian matrix in this case is a triangular matrix. As a result, the Jacobian determinant is the product of the diagonal entries:

The inequality in Eq 8 is thus:

where $\mathbf{Y}=\{Y,X_{2},X_{3},\dots,X_{d}\}$.

On the left hand side of the above inequality, the joint entropy of $\mathbf{Y}$ can be expressed using the conditional entropies as $H(\mathbf{Y})=H(Y_{1}|\mathbf{X}_{2\sim d})+H(\mathbf{X}_{2\sim d})$ [21]. The subscript $2\sim d$ indicates the index ranging from $2$ to $d$. On the right hand side of Eq 10, we have $H(\mathbf{X})\leq H(X_{1})+H(\mathbf{X}_{2\sim d})$ using the subadditivity property of the joint entropy of the input variables. The joint entropy $H(\mathbf{X})$ becomes additive if the input variables are independent.

Putting these together, the inequality of Eq 10 then becomes:

where the expectation is with respect to the jPDF of the input variables. Note that $H(Y|\mathbf{X}_{2\sim d})=\mathbb{E}\left[H(Y|\mathbf{x}_{2\sim d})\right]$ is the expected conditional entropy.

The reasoning above uses the first variable $x_{1}$ as an example. However, the results hold for any variables via simple row/column exchanges, which only affects the sign of the determinant but not its modulus.

Therefore, we have the following equality:

where $\sim i$ indicates the index ranges from $1$ to $d$ excluding $i$, and we have introduced the following notation $l_{i}=\mathbb{E}\left[\ln{\left|{\partial g(\mathbf{x})}/{\partial x_{i}}\right|}\right]$.

The conditional entropy $H(Y|\mathbf{X}_{\sim i})$ measures the remaining entropy of $Y$ if the true values of $\mathbf{X}_{\sim i}$ can be determined, in analogy to the total effect index $S_{T_{i}}$ from variance-based indices. Note that the conditional entropy $H(Y|\mathbf{X}_{\sim i})=\mathbb{E}\left[H(Y|\mathbf{x}_{\sim i})\right]$ is a global measure of uncertainty as the expectation is with respect to all possible values of $\mathbf{X}_{\sim i}$.

The above inequality thus demonstrates that, for a differentiable function $y=g(\mathbf{x})$, the entropy-based total sensitivity is bounded by the expectation of log partial derivatives of the function, with the addition of the entropy of the input variable of interest. And this inequality becomes an equality if the input variables are independent and the transformation $g(\cdot)$ is bijection, i.e. $g:X\rightarrow Y$ is an invertible transformation and it has a unique inverse.

The inequality in Eq 11 is one of the main contributions of this paper. It establishes an upper bound for the total effect entropy $H(Y|\mathbf{X}_{\sim i})$ using computationally efficient partial-derivative based functionals. As smaller $l_{i}+H(X_{i})$ tends to indicate smaller total effect entropy, it can thus be used to screen uninfluential variables and work as a low cost proxy for entropy-based indices.

The use of $l_{i}+H(X_{i})$ as a proxy for the total effect entropy is similar to the DGSM-based upper bound for the variance-based $S_{T_{i}}$ described in Eq 5. However, there is a major issue with the inequality in Eq 11: we cannot normalise Eq 11 by $H(Y)$ to obtain a bound on the entropy-based total effect indices $\eta_{T_{i}}$, as $H(Y)$ can potentially be negative. In Section 4, we will overcome this issue by taking the exponentiation of the inequality in Eq 11, and demonstrate that the exponential entropy based measure has a more intuitive sensitivity interpretation and is closely linked to the variance-DGSM proxy. Before that, we will first verify Eq 11 using several examples in the next two subsections.

We will first test monotonic functions where the equality of Eq 11 is expected, and examine the inequality with more general functions in the next section. All input variables are assumed to have the same uniform distribution for examples 1 - 4, i.e. $x_{i}\sim\mathbb{U}(0,1)$, while Gaussian distributions are used for example 5. For verification purposes, all the examples in this section are chosen to have tractable expressions for both the integral of derivatives and the conditional entropies. Note that for simplicity the input variables are assumed to be independent in these numerical examples. Although input independence is required for the equality condition, the inequality in Eq 11 makes no assumptions of independence or distribution type for $X_{i}$.

Both Ishigami function and G-function are commonly used test functions for global sensitivity analysis, due to the presence of strong interactions. These two functions, each with three input variables, are used in this section to demonstrate the inequality relationship derived in Eq 11. We will focus on the verification of Eq 11 in this section, but will re-visit these two functions in the next section for a discussion in comparison with variance-based sensitivity indices.

The conditional entropies are estimated numerically using Monte Carlo sampling as described in Appendix A. Different numbers of samples are used, ranging from $1e6$ to $1e8$. For each estimation, the computation is repeated 20 times and both the mean value and the standard deviation (std) are reported in Table 4 and 5. For the estimation of derivative-based $l_{i}$, Matlab’s inbuilt numerical integrator "integral" is used with the default tolerance setting.

This becomes more evident if we plot the histogram of function output $Y$ from examples 1 to 3 in Figure 4. These histograms are based on $10^{6}$ Monte Carlo samples and from which we have calculated the exponential entropy $e^{H(Y)}$ for each case. When compared to the range of the distribution, the values of the exponential entropy provide intuitive indication of the spread of the distribution, in analogy to standard deviation or variance. For example, the average width of the top and bottom of the trapezium in Example 3 is approximately $(3+4)/2=3.5$, which is about the same as the exponential entropy of $3.543$.

Therefore, the exponential entropy $e^{H}$ can be regarded as a measure of the extent, or effective support, of a distribution and this has been discussed in details in [6]. As $H(Y|\mathbf{X}_{\sim i})$ measures the remaining entropy in average if the true values of $\mathbf{X}_{\sim i}$ can be determined, the exponentiation of the conditional entropy, $e^{H(Y|\mathbf{X}_{\sim i})}$, can thus be regarded as the effective remaining range of the output distribution conditioning on that $\mathbf{X}_{\sim i}$ are known. $\kappa_{T_{i}}$ then measures the ratio of the effective range before and after $\mathbf{X}_{\sim i}$ are fixed, and larger $\kappa_{T_{i}}$ thus indicate a higher influence of $X_{i}$.

Last but not least, in addition to its non-negativity and a more intuitive interpretation for GSA, exponential entropy is also closely linked to variance-based GSA indices and their corresponding bounds. To demonstrate this, we first note that the three different derivative-based sensitivity indices listed in Table 2 are closely related as:

where it is evident that $\mu_{i}\leq\sqrt{\nu_{i}}$ based on Cauchy-Schwarz inequality. In addition, we have $e^{l_{i}}\leq\mu_{i}$ using Jensen’s inequality as the exponential function is convex. So the inequality for the exponential entropy in Eq 12 can be further extended to:

where we recall that $\nu_{i}$ is the derivative-based DGSM indices.

Eq 15 already looks remarkably similar to the variance-DGSM inequality given in Eq 5. In fact, the squared exponential entropy $e^{2H(X)}$ of a random variable $X$, also called entropy power, is known from information theory to be bounded by the variance of $X$ with $e^{2H(X)}\leq 2\pi e{V(X)}$ [21], where the equality is obtained if $X$ is a Gaussian random variable. Therefore, for independent Gaussian inputs with variance $\sigma_{i}^{2}$, we have $e^{2H(X_{i})}=2\pi e\sigma_{i}^{2}$. The squared exponential entropy $e^{2H(X)}$ thus measures the ‘effective variance’ as it is simply the variance of a Gaussian random variable with the same entropy [25].

In the special case where the function $y=g(\mathbf{x})$ is a linear function, the conditional distribution of $Y|\mathbf{X}_{\sim i}$ is also Gaussian, and it indicates that the conditional entropy power would also be proportional to the conditional variances of the output, i.e. $e^{2H(Y|\mathbf{X}_{\sim i})}=2\pi eV(Y|\mathbf{X}_{\sim i})$ (cf Example 5 in Section 3). Therefore, it is clear that with independent Gaussian inputs and a linear function, we have from Eq 15:

which indicates that in this special case the entropy-DGSM bound from Eq 15 is equivalent to the variance-DGSM relationship given in Eq 5.

In comparison, we have also calculated the variance based $S_{T_{i}}$, using $10^{5}$ samples with 20 repetitions. The mean values for $S_{T_{i}}$ are listed in Table 6. It can be seen that the results from $\kappa_{i}$ are generally consistent with $S_{T_{i}}$, especially for G-function where the sensitivity ranking are similar both qualitatively and quantitatively. For the Ishigami function, both indices have successfully identified the contribution of $X_{3}$ which is the lowest. However, the relative importance of $X_{1}$ and $X_{2}$ are opposite from $\kappa_{i}$ and $S_{T_{i}}$. This difference between variance-based and entropy-based ranking for Ishigami function was also noted in [3], but without further explanation.

We attempt to look at the difference using a simple function with $(x,z)$ as the input random variables:

which can be seen to represent the interaction effect of $x_{1}$ and $x_{3}$ in the Ishigami function. The conditional variance and entropy of $Y|Z=z$ are thus:

where the translation invariance of variance and entropy are used. And the corresponding un-conditional $V(Y|Z)$ and $H(Y|Z)$ can be found as:

From 19, it is clear that despite the similarities, the conditional variance $V(Y|Z)$ depends on $\alpha^{2}$ which grows much faster than $\ln|\alpha|$ for the conditional entropy $H(Y|Z)$ as $|\alpha|$ increases.

In the Ishigami function, the effect of $x_{3}$ grows fast as $x_{3}$ approach its boundary, i.e. towards $-\pi$ and $\pi$. This is evident from the scatter plot of for $x_{3}$ in Figure 3a. The much stronger interaction at the support boundary between $x_{1}$ and $x_{3}$ has also been noted in [26], where the importance of $x_{1}$ was found to be much reduced if the support of the distribution of $x_{3}$ is reduced by $10$ percent. Putting $x_{3}\sim[-\pi+\pi/10,\pi-\pi/10]$ would put $x_{2}$ as the most influential variable.

From Eq 19, we can see that the interaction between $x_{1}$ and $x_{3}$ is more influential for the conditional variance due to the squaring effect, as compared to the entropy operation which takes logarithm of the interaction. This difference increases towards the boundary as the interaction between $x_{1}$ and $x_{3}$ gets stronger towards $-\pi$ and $\pi$. And this helps to explain why $x_{1}$ is the most influential variable for the variance-based $S_{T_{i}}$. This example highlights that, despite many similarities, entropy and variance are fundamentally different, for example, the variable interactions are processed differently between them.