SOBOL' INDICES FOR PROBLEMS DEFINED IN NON-RECTANGULAR DOMAINS

. A novel theoretical and numerical framework for the estimation of Sobol’ sensitivity indices for models in which inputs are confined to a non-rectangular domain (e.g., in presence of inequality constraints) is developed. Two numerical methods, namely the quadrature integration method which may be very efficient for problems of low and medium dimensionality and the MC/QMC estimators based on the acceptance-rejection sampling method are proposed for the numerical estimation of Sobol’ sensitivity indices. Several model test functions with constraints are considered for which analytical solutions for Sobol’ sensitivity indices were found. These solutions were used as benchmarks for verifying numerical estimates. The method is shown to be general and efficient.


Introduction
High complexity of models in physics, chemistry, environmental studies, reliability and structural analysis and other fields results in the increased uncertainty in model parameters and model structures. Uncertainty in inputs is reflected in uncertainty of model outputs and predictions. Uncertainty and sensitivity analysis has been recognized as an essential part of model applications. Global sensitivity analysis (GSA) offers a comprehensive approach to model analysis by quantifying how the uncertainty in model output is apportioned to the uncertainty in model inputs [1,2]. Unlike local sensitivity analysis, GSA estimates the effect of varying a given input (or set of inputs) while all other inputs are varied as well, thus providing a measure of interactions among variables. GSA is used to identify key parameters whose uncertainty most affects the output. This information then can be used to rank variables, fix unessential variables and thus decrease problem dimensionality. The variance-based method based on Sobol' sensitivity indices has become very popular among practitioners due to its efficiency and ease of interpretation [3,4]. Most of the developed techniques for GSA are designed under the hypothesis that inputs are independent. However, in many cases there are dependencies among inputs, which may have significant impact on the importance results.
There have been a number of attempts to extend GSA to models with dependent inputs. We present a brief survey of only the most recent and easy-to-follow developments. Xu and Gertner [5] suggested to split the contribution of an individual input to the uncertainty of the model output into two components: the correlated contribution and the uncorrelated one. A regression-based method was used for estimating the correlated and uncorrelated contributions of the inputs. The approach developed originally only for linear models was later extended to nonlinear models [6]. Li et al [7] proposed a generalization of the ANOVA-HDMR decomposition by including covariances to the model variance. They distinguished between structural and correlative contributions of a given input. This method presents some critical issues such as non-uniqueness of the functional decomposition and hence difficulties in interpreting the results. Da Veiga et al. [8] used the estimator of the main effect index firstly proposed in [9] and further developed in [10]. They suggested a new procedure based on a heteroscedastic regression model and local polynomial metamodel.
A novel approach for the estimation of Sobol' sensitivity indices for models with dependent variables using generalizations of direct Sobol' formulas was developed in Kucherenko et al. [11]. Both the main effect and total sensitivity indices were derived as generalizations of Sobol' sensitivity indices. Formulas and Monte Carlo (MC) numerical estimators similar to Sobol' formulas were proposed. A Gaussian copula-based approach was used for sampling from arbitrary multivariate probability distributions. The generalization does not involve the use of surrogate models, data-fitting procedures or orthogonalization of the input factor space.
Mara and Tarantola [12] introduced a set of sensitivity indices which relies on orthogonalisation of correlated inputs. The computation of sensitivity indices was performed using a parametric method, specifically the polynomial chaos expansion. Mara et al. [13] extended the development of sensitivity indices suggested in [12] by proposing two sampling strategies for non-parametric, numerical estimation using the Rosenblatt transformation (RT) [14]. RT, and hence values of sensitivity indices, is not unique: for a model with n inputs there are n! possibilities corresponding to all possible permutations of the elements of the input vector 1 ,..., n x x x . The authors considered only the RT obtained after circularly reordering the set 1 ,..., n xx , resulting in n RT transformations. The authors also established the link with the indices proposed by Kucherenko et al. [11]. W.Hao et al. [15] suggested a detailed interpretation of indices proposed in [11]. Taking a quadratic polynomial model with normal inputs as an illustration, they derived explicit expressions for sensitivity indices and considered contributions to the values of sensitivity indices from all components. In this work we propose GSA formulations for an even wider class of problems with dependent variables, namely for models involving inequality constraints (which naturally leads to the term 'constrained GSA' or cGSA). Such constraints impose structural dependences between model variables in addition to potential correlations between them. This implies that the parameter space may no longer be considered as an n-dimensional hypercube or hyperrectangle as considered within GSA so far, but may assume any shapes (including disconnected domains) depending on the number and nature of constraints. This class of problems encompasses a wide range of situations encountered in the natural sciences, engineering, design, economics and finances where model variables are subject to certain limitations imposed e.g. by conservation laws, geometry, costs, quality constraints etc. An important particular case within this paradigm corresponds to imposing a minimum (maximum) threshold for the model output, i.e., The development of efficient computational approaches for cGSA is challenging because of potentially arbitrary shape of the feasible domain of the model variables' variation, thus requiring the development of special MC or quasi MC (QMC) sampling techniques and approaches for computing sensitivity indices. This becomes especially difficult for models of high dimensionality. Another challenge consists in analysing and interpreting model variance and sensitivity indices in such domains, since in this case sensitivity indices will carry structural information imposed by model constraints. The interpretation of sensitivity indices in such circumstances is crucial for the efficient design of experiments as well as for potential model reduction by fixing or eliminating nonessential variables.
In this paper we have developed and tested several approaches for the numerical estimation of main effect and total sensitivity indices in the cGSA setting. It is organized as follows: The next Section presents formulas for the main effect and total sensitivity indices for models with dependent variables and acceptance-rejection method which can be used to avoid sampling from conditional distributions. Section 3 considers numerical aspects of the method, namely it presents general MC estimators, acceptancerejection estimation of Sobol' sensitivity indices using grid quadrature formulas and MC estimators. Section 4 presents the application of the proposed methodology to three test cases. Numerical and analytical results are compared and the convergence rates of the MC and QMC methods are discussed. Finally, conclusions are presented in the last Section.

Variance-based sensitivity indices for models with dependent inputs
The ratio is known as the main effect index of the subset y . The quantity is known as the total effect index of the subset y . Collectively y S , T y S are known as Sobol' indices. In this paper our notations follow closely those of Kucherenko et al. [11].
( ( , )) z E f y z in Eq. (1) denotes a conditional expectation with respect to z with y being fixed and ( ( , )) y D f y z in (2) is a conditional variance with respect to y with z being fixed. We use notations z and z (correspondingly y and y ) to distinguish a random vector z (correspondingly y ) generated from a joint pdf p y z and a random vector z (correspondingly y ) generated from a conditional pdf ( , | ) p y z y (correspondingly ( , | ) p y z z ).
A formula for the main effect index can be explicitly written as ( Here py is the marginal pdf of subset of inputs y . This is the so-called double loop formula (DL) [1]. It can be transformed into a different formula which in some cases can be more efficient [11]: (4) The following notation is used: z , z are two different random vectors generated from the joint pdf ( , ) p y z , a random vector ẑ is generated from a conditional pdf ( , | ) p y z y .
A formula for computing the total effect index T y S (2) can be written as

Transformation from conditional to joint and marginal distributions
We assume that there is a known procedure for the generation of vectors ( , ) yz from a joint pdf ( , ) p y z defined in a non-rectangular domain n , which means in particular that ( , ) 0 p y z Here the lower-dimensional integrals are understood in the sense that The usefulness of these formulas is based on the ability to sample from the joint pdf ( , ) p y z in n . Although there is a technique based on the sequential sampling from an inverse cumulative distribution presented in Appendix A, it has a limited applicability. One example of an application of this technique is given in Appendix B. In the next section we present a more general approach based on the acceptancerejection method.

Acceptance -rejection method
Consider the joint pdf ) , ( z y p in the absence of constraints in n H . We assume that constraining the variables to an area nn H implies that their joint pdf becomes n n yz I y z yz (12) The latter can be represented through the Heaviside step function () U applied to the constraint(s): ( , ) ( ( , )) I y z U g y z .
A constrained marginal pdf is then defined as correspondingly.
Using the expressions above for pdf's as well as the expected value and total variance of the function, the main effect and total indices can be computed through their integral formulations presented in (7)-(9). The important difference, however, is that even if the relevant pdf's are known only within the enveloping hypercube (i.e., for the unconstrained case), they can be used to directly compute their constrained counterparts on the basis of an acceptance-rejection approach invoking the indicator function of the feasible subdomain.
The DL formula (7) can be rewritten in the following simplified form An alternative expression for the total effect index T y S can be obtained using the identity Using the Bayes formula and transformations similar to presented above we obtain a DL-like formula for total indices:

Interpretation of indices
We start with Interpretation of the total index. In the case of independent inputs where yz D is an interaction term between subsets y and z . From definition (19) it follows that T y D is the part of D which remains after deduction of ( , ) zy D E f y z . In the case of dependent inputs zy D E f y z contains contribution to the total variance corresponding to the effect of subset z on its own and its dependence with subset y due to correlation or dependence via inequality constrains, which we'll call for brevity dependence contribution

Numerical methods
In this section we present general MC estimators, acceptance-rejection estimation of Sobol' sensitivity indices using grid quadrature formulas and MC estimators based on the acceptance-rejection method.

MC estimators
We consider MC estimators for the evaluation of integral expressions in (7)- (9) assuming that all pdf's are defined in n and there is a sampling procedure able to generate random vectors within this domain. Formula (7)  (23) The subdivision into bins is done in the same way for all inputs using the same set of sample points. A critical issue is the link between N and y N . It was suggested in [16] to use as a "rule of thumb" . y NN We note that it is practically impossible to generalise the DLR formula for more than one index, hence there is no similarly efficient DLR formula which allows to compute total Sobol' sensitivity indices.
The expected value and total variance are computed using the MC estimators and an estimator for the total effect index based on formula (9) can be written as: The usefulness of these estimators relies on the ability to sample from the joint pdf ( , ) p y z in n .
In the next two sub-sections we present a practical approach for such sampling based on the acceptancerejection method.

Acceptance-rejection estimation of sensitivity indices using grid quadrature formulas
Integrals in (16)-(18) and (21) can be efficiently evaluated using grid quadrature methods for problems of low and medium dimensionality. Application of grid quadrature methods is rather straightforward when integration domains are hyperrectangular as in the case of the acceptance-rejection approach presented in section 2.3. However, given that the feasible domain n can be arbitrary in shape, efficient higher-order quadrature methods (such as those based on Gaussian or Clenshaw-Curtis quadrature [17]) lose their advantages as the integrands are not differentiable at the boundary of the permissible domain n (the indicator function has a jump discontinuity at the boundary of n ). Thus the use of lower-order integration formulas such as the second-order multidimensional trapezoidal rule used in this paper is preferable because they are less sensitive to non-smooth or discontinuous integrands.
Numerical integration can gain additional efficiency by performing preliminary bracketing of the domain n (or finding its 'minimum bounding box') within n H to minimize the number of rejected points during sampling. This is especially important in higher dimensions, when the number of grid points in each dimension can't be large due to the "curse of dimensionality". The total number of grid points is n k N , where k is the number of points in each dimension. k includes both internal and boundary points and it can't be smaller than 4. We don't consider the case of the different values of k for different directions due the lack of prior information about importance of individual inputs in the general case.
where i x is the vector of all model variables but This approach was successfully implemented for problems of dimensionality up to 10. The results are presented in Section 4. The use of this approach for higher dimensional models is computationally prohibitive. This is also the main reason why the grid quadrature approach is not applicable to the modified formulas (8)-(9) since the effective number of dimensions for integration in this case is

MC estimators for the acceptance-rejection method
In this subsection we consider MC estimators of y S and T y S valid when ( , ) p y z is known but the constrained pdf ( , ) p y z is not known explicitly. We also assume that there is no explicit technique for sampling in n of arbitrary shape. Firstly we derive an MC estimator for y S defined by (18). Its approximation leads to the DLR approach discussed in subsection 3.1. Based on a sequence of points i.e., by averaging within each of y N bins, where j is the index of a bin and A j y represents an average value of y in the j -th bin. The marginal distribution function () py in the denominator of (18) is approximated within each bin by Combining the above we obtain a DLR MC estimator for the main effect index defined by formula (18) in the following form: An MC estimator for y S based on the modified formula (8) has the form: Similarly an estimator for the total effect index (9) can be written as:  (38) 11 We note that only the indicator function is evaluated at each of these points and although this additional sampling at each l z increases the computational cost of the estimator (37), this increase may not be detrimental if the computational cost of the evaluation of model constraints is significantly smaller than that of the model function itself (which is often the case in practical problems). A value 64 z N was used in the test cases reported below, which was found to be high enough to provide good convergence rates.

Test cases
For all test cases considered in this section we found analytical values of sensitivity indices, so that they can be used as benchmarks for verification of numerical estimates.

Product function in triangle domain
1 .

Consider the function
Expectation and total variance for this function can also be computed explicitly: (44) Using the definition (7) 80  1  5  7  2  3 12 27

g-function in triangular domains
Consider the so-called g-function which is often used in GSA for illustration purposes [1]: with parameters 0 1 a and 1 2 a (see Fig. 3) and 12 , xx uniformly distributed in the unit square. Below we compute its sensitivity indices when the feasible domain is restricted by a parametric linear (Fig. 4) or nonlinear (Fig. 9) constraint.

A linear constraint
In this subsection we consider a more general parametric linear constraint (illustrated in Fig. 4) than that used in the test case presented in the previous subsection. The shape and size of the feasible domain 2 are defined by the variable angle between the top side of the unit square and the line defined by the linear constraint Firstly we obtain analytical results for two values of : 6 and 4 . While the latter value leads to a symmetrical domain (complementary to 1 in Fig. 2), the former does not result in any specific simplification of the problem. The reference values of the mean, total variance and main effects are given in Table 1. Exact solution for 6 was obtained using symbolic integration in Maple®, however the resulting expressions are too cumbersome to report them here, so the values are given in the decimal form to sufficient precision for error analysis (see below). As noted above, the values of total sensitivity indices are readily obtained from the relationships  Table 1 are denoted by overlaid symbols. The results perfectly respect limiting cases given in Table 1. The numerical results were obtained with the use of the grid quadrature method.    We also tested the performance of the MC estimators. The root-mean-square error (RMSE) obtained using the MC estimator (35) for the main effect indices and the DLR estimator (37) for total effects indices for a fixed value of 6 is presented in Fig. 7. To reduce the scatter in the error estimation the values of RMSE were averaged over L = 50 independent runs:  Finally, we compare the performance of the two presented estimators for the main effect indices: namely, the DLR (35) and the estimator (36) for the modified formula (8) by S. Kucherenko et al [11], the latter is denoted as SK. The results presented in Fig. 8 show significant advantages of the DLR approach when considering the convergence rate versus the total number of sampling points CPU N .

A parabolic constraint
In this subsection we consider a more complex example involving a nonlinear constraint of the form: which describes the part of the unit square above a parabola (Fig. 9). Depending on the value of the permissible domain is either connected ( 4 ) or disconnected ( 4 ) as illustrated in Fig. 9. We note that the domain is non-convex for any value of .

K-function
K-function is defined as These constraints can be represented using the following indicator functions: (1 ) (1 ) I U x x . See Fig. 13 for a schematic plot illustrating 1 I constraint in the 3D space.
For numerical estimates we make use of grid quadrature formulas (multidimensional trapezoidal rule) presented in subsection 3.2. In order to assess the accuracy of numerical computations in the unconstrained case the exact solution for the total effect indices reported in [19] was used while the analytical solution for the main effect sensitivity indices was derived in this paper:

Conclusions
In this work, we have proposed a novel concept of constrained GSA which adds the ability to analyse model output variance in arbitrarily shaped n-dimensional domains. This amounts to greatly expanding the scope of GSA by allowing model variables to be subject to inequality constraints, which is common in a range of situations of practical importance.
The proposed formulas build upon Sobol' sensitivity indices and their recent development for models with dependent variables [11]. The advantage of the presented formulations is that no prior knowledge of conditional or marginal distributions is assumed. All the required dependences are derived from the joint pdf in the presence of constraints. It is shown that the knowledge of the joint pdf corresponding to the unconstrained formulation is sufficient to build numerical estimators for sensitivity indices.

22
Three types of numerical estimators of the sensitivity indices have been proposed. Grid quadrature may compete with MC estimators for low and medium dimensional models. However, its convergence rate rapidly degrades with increasing model dimensionality. Despite its simple concept, the DLR approach demonstrates good convergence when applied to the evaluation of main effects sensitivity indices outperforming the MC estimator of the modified Sobol' formulas. On the other hand, DLR is not a viable approach for the evaluation of total effect indices for which the modified Sobol' formula gives good results.
Further work is needed to develop a clear interpretation of the cGSA results, in particular by decomposing the variance contribution into correlated and structural (uncorrelated) contributions.