# Asymptotic Analysis for the Variance-Based Global Sensitivity Indices

Authors: Pavel M. Bokov
Journal: Science and Technology of Nuclear Installations
 (2012)

---

## Abstract

We discuss the estimation of the uncertainty and sensitivity parameters for a model response under the assumption that the input variables are normally distributed and block-wise correlated with the covariance matrix, which is small in some norm. These conditions may arise when considering the impact of the group-wise neutron cross-sections' uncertainties on the uncertainty of some reactor parameters such as the neutron multiplication factor. The variance-based global sensitivity analysis, considered in our work, involves the calculation of multidimensional integrals. When the input uncertainties are small, the values of these integrals can be estimated using an asymptotic analysis method called the Laplace approximation. The asymptotic formulas for the output variance and for the global sensitivity indices have been obtained using the Laplace approximation method. It is demonstrated that the asymptotic formula for uncertainty propagation matches the uncertainty propagation formula being used in the local sensitivity analysis. The applicability of the obtained asymptotic approximations was successfully demonstrated on a test problem with realistic cross-section and covariance matrix values.


---

## Body

## 1. Introduction

Uncertainty analysis (UA) is a mathematical tool that allows one to quantify the uncertainty of the model output as a result of the uncertainty in the model input. If the input consists of more than one variable, another mathematical tool, sensitivity analysis (SA), may be used to quantify the contribution of each parameter to the output uncertainty. SA and UA methods can be gathered in two broad families: local sensitivity analysis (LSA) and global sensitivity analysis (GSA).LSA methods allow one to analyze the behavior of the model output in the vicinity of a chosen point. They are usually efficient in computer time but may be inadequate for nonlinear models. GSA methods allow one to explore the full-phase space of input parameters and to take the nonlinearity of the model into account. GSA methods are generally more computationally intensive than LSA methods. The major difference between these two approaches can be summarized with the following quote from Sobol’ [1]: “Global sensitivity indices should be regarded as a tool for studying the mathematical model rather than its specified solution.” A detailed discussion on SA methods can be found in [1–3]. The overview of SA as applied in the nuclear reactor calculation field is given in [4] and the references therein.In this paper we will concentrate our attention on the so-called variance-based GSA. Amongst different GSA methods, the variance-based methods are the methods of preference because they are well established, robust and model-independent. On the other hand, despite the attractiveness of the variance-based GSA, UA, and SA for nuclear reactor problems are traditionally based on the LSA methods and rarely on GSA methods. The reason why the variance-based GSA method was not widely applied is that there are a few obstacles in the way of their implementation in reactor calculations. The first obstacle is the correlation between input parameters, that makes the variance-based GSA in its traditional formulation meaningless. This limitation may appear, for example, when considering multigroup neutron cross-sections as input parameters, which may have correlations between cross-sections from different energy groups. The partial remedy for this problem was proposed by Jacques et al. [5], who remarked that the variance-based GSA is still applicable to the mutually independent subsets of input variables. Another, more radical solution, is to abandon the variance-based GSA and to use, for instance, entropy-based GSA [6, 7]. Entropy-based GSA methods can be used for correlated inputs, but they are not as well-established as variance-based ones, are algorithmically more complex and computationally more expensive. The second obstacle is the potential impracticality (and even intractability) when the number of input parameters is big. This problem appears because GSA methods involve the calculation of multidimensional integrals and is traditionally circumvented by applying the integration methods suitable for multidimensional problems, such as Monte Carlo [1, 8] or sparse grid [9, 10] quadratures. These quadratures are often applied in conjunction with model reduction methods [11–13], which are often based on the functional ANOVA (analysis of variance) decomposition [8, 14]. When applied in the field of nuclear reactor calculations, UA and SA methods have an additional limitation: the model value is often the result of an expensive calculation and this limits the number of samples available for the analysis.Nevertheless, there are some indications that a large number of samples may not be needed for practical applications. Adetula and Bokov [15, 16] described a numerical method for the calculation of global sensitivity parameters for multigroup cross-section dependent problems with block-wise correlated inputs. The method was illustrated using a test two-group problem that involved realistic cross-section data. Various integration techniques were used and one of them, sparse grid quadrature, has demonstrated a surprisingly good accuracy with a small number of sample points. Moreover, it was observed that the interaction between inputs is negligible. Hence, it was concluded that the model is virtually additive and can be approximated with low-order polynomials. This result was explained with the following hypothesis: the uncertainty in cross-sections is so small that only the vicinity of the cross-section mean values contributes to the integrals used in the estimation of sensitivity indices. In this vicinity the neutron multiplication factor can be approximated with a fairly linear function. This hypothesis was tested and confirmed empirically (see [15, 16] for detail), but none theoretical support has been provided.The present work aims to fill this gap and to provide a theoretical explanation and a quantitative description of this phenomenon. The idea behind this hypothesis is rather simple and may be illustrated with the sketch given in Figure1 for a one-dimensional case. Consider a model represented with a smooth function f(x), where the input variable x is normally distributed with a mean μ and a standard deviation σ. Let us assume that it is necessary to evaluate the integral in the following form:
(1)I=∫-∞+∞‍f(x)p(x)dx,
where p(x) is the probability density function for x. As one can observe in Figure 1, only the interval under the Gaussian (which is of the order of a few σ) effectively contributes to the integral. In the limit, when σ→0, that corresponds to the case σ=σ0, the probability distribution becomes the delta function and the contribution of the function to the integral is point-wise, that is, I=f(μ). When σ is small enough, as in the case σ=σ1, the part of the function f(x), contributing to the integral, is almost linear. Finally, if the the Gaussian spans wider intervals, then the nonlinearity of f(x) becomes progressively more prominent (see cases for σ=σ2 and σ=σ3). A similar behavior can be observed in the multidimensional case. The qualitative description of this behavior in a one- and multidimensional case can be given using a method known as the Laplace approximation through the asymptotic approximation of the involved integral.Figure 1
Effective vicinities of the cross-section mean, which contribute to the integral, for the standard deviationsσ3>σ2>σ1>σ0=0. Dots on the curve y=f(x) represent the boundaries of the corresponding effective intervals on which integrand is essentially not zero.In this work we intend to apply the Laplace method for the calculation of the output variance and sensitivity indices as they are defined in GSA. The application of the Laplace approximation in a one-dimensional case is relatively straightforward, while in a multidimensional case it may become cumbersome, this is why we limit our approximation with the first nonvanishing term in the approximation. Nevertheless, as we will demonstrate below, this approximation is sufficient to explain the results as obtained in [15, 16].The rest of the paper is organized in the following way. Section2 provides necessary definitions followed by a derivation of the asymptotic approximation. Section 3 is used to discuss the asymptotic approximation. Section 4 contains the description and analysis of the test problem. Finally, Section 5 is used to present our conclusions.


## 2. Asymptotic Approximation of Sensitivity Indices

### 2.1. Variance-Based Global Sensitivity Analysis

Let𝒳={X1,X2,…,Xd} be a set of continuous random variables, 𝒳u be its arbitrary subset and 𝒳v be a subset of complimentary variables, that is, 𝒳=𝒳u∪𝒳v. The random variables can be gathered into the corresponding column vectors: Xu, Xv, and X=(Xu,Xv). Following the arguments given in [5, 15, 16], we assume that random variables from 𝒳u and 𝒳v are mutually statistically independent, hence their joint probability function is the product of the marginal probability density functions:
(2)p(x)=p(xu)p(xv).
Note that we will use, as is the rule in statistics, a capital letter (e.g., Xi) to denote a random variable and the corresponding lowercase letter (xi) to denote its value (realizations). Let us assume that the random variables Xi(i=1,2,…,d) are distributed according to the normal law with known means and covariances. The multivariate normal distribution is characterized by the probability density function [17]
(3)p(x)=1(2π)d/2det(Σ)1/2exp[-12(x-μ)TΣ-1(x-μ)],
where X is the column vector of the input random variables, μ=E[X] is the column vector of their expected values, Σ=E[(x-μ)(x-μ)T] is the covariance matrix and E denotes the mathematical expectation and the symbol “T” denotes the operation of transposing a row to a column.Consider a model(4)Y=f(X1,X2,…,Xd)=f(X),
where the function f:ℝd→ℝ is generally nonlinear and may result from numerical calculations. Xi are called inputs and Y is called the output or response. The response Y=f(X) is a random variable itself, with the expected value f∅=E[f(X)] and the variance D=E[(f(X)-f∅)2]. The variance of the output, D, can be used to characterize the uncertainty in the response due to the uncertainty in the input, or in other words, the uncertainty propagation. Note that D can be calculated, regardless of whether the normally distributed inputs are correlated or not. Sensitivity indices characterize the contribution of a subset 𝒳u (and, correspondingly, vector Xu) to the response variance and they are meaningful only if subset 𝒳u and its complementary subset 𝒳v are statistically independent. Two types of global sensitivity indices, with respect to a subset 𝒳u, can be introduced:(1)
SXu, which is called the main sensitivity index (MSI) and it represents the effect due to Xu only;(2)
SXutot, which is called the total sensitivity index (TSI) and it represents the contribution to the variance of Xu along with all the interactions of this variable with other variables.The global sensitivity indices can be calculated using the ratios [1, 8]
(5)SXutot=DXutotD,SXu=DXuD,
where parameters D, DXu, and DXutot can be calculated through the following multidimensional integrals [12, 15, 16]:
(6)f∅=∫ℝd‍f(x)p(x)dx,(7)D=∫ℝd‍f2(x)p(x)dx-f∅2,(8)DXu=∫ℝ2d-dim(Xu)‍f(x)f(xu,xv′)p(x)p(xv′)dxdxv′-f∅2,(9)DXutot=12∫ℝd+dim(Xu)‍[f(x)-f(xu′,xv)]2p(x)p(xu′)dxdxu′,
which are a generalization of the classical formulas given by Sobol’ [1, 8]. Here the prime symbol over a variable (e.g., as in xu′) means that this variable has to be sampled independently from the corresponding marginal distribution (p(xu′) in this case) of its unprimed analogue and dim(·) is the number of elements in a vector, that is, the dimensionality of the vector. Calculation of the sensitivity indices requires the calculation of integrals (6)–(9), which can be written in the following general form:
(10)E[g]=∫ℝdeff‍g(x~)p(x~)dx~,
where the integration operator E[·] denotes the mathematical expectation with respect to extended random variables X~, g(x~) represents a function being integrated and deff=dim(X~) is the effective dimensionality of the integral. For instance, in integral (8) function g(x~) represents [f(x)f(xu,xv′)], x~ represents (xu,xv,xv′), and the effective dimensionality is deff=2d-dim(Xu). Similarly, in the integral (6) function g(x~) represents f(x), x~=x and deff=d.


### 2.2. Laplace Approximation of Integrals

In the limit when elements of the covariance matrix are small in some norm, one can use for the calculation of integrals (6)–(9) the asymptotic analysis technique known as the Laplace approximation. The Laplace approximation is an approximation to integrals in the following form:
(11)I(λ)=∫Ω‍h(x)exp[λR(x)]dx,
as λ→∞ with λ>0. Here Ω is a measurable subset of ℝd and R(x) is a real-valued function of x. When applying the Laplace approximation formula, it is assumed that R(x) is an infinity differentiable function with its maximum value at a single point x0 in the interior of Ω (see [18] for details) and h(x) is a sufficiently smooth function of its argument.Integrals in the form (10), used for the calculation of sensitivity indices, fulfill the above requirements for the applicability of the Laplace approximation formula if one considers x0=μ and assumes that elements of the covariance matrix Σ are small. Let us assume that the covariance matrix Σ is small in some norm, that is, in the limit ∥Σ∥→0. For the sake of the asymptotic analysis, let us introduce the big parameter λ as the inverse absolute value of the determinant of the covariance matrix
(12)λ=|det(Σ)|-1.
We will assume that the covariance matrix is strictly positive definite, hence the absolute value sign in (12) may be omitted. Factoring out the coefficient λ allows us to introduce the normalized covariance matrix
(13)C=Σ|det(Σ)|-1=λΣ.
Taking into account that det(C)=det(λΣ)=λddet(Σ), the probability density function (3) can be written in a form suitable for the Laplace approximation:
(14)p(x)=1|det(C)|1/2(λ2π)d/2exp[-λ12(x-μ)TC-1(x-μ)].
The covariance matrix Σ is real-valued, symmetric, and positively definite, hence matrix C allows the Cholesky factorization:
(15)C=LLT,
where L is a lower triangular matrix with strictly positive diagonal entries. The exponent in the joint probability density function can then be transformed in the following way:
(16)λ12(x-μ)TC-1(x-μ)=λ12zTz,
where new random variables Zi(i=1,…,d) are introduced by the linear transformation z=L-1(x-μ). The initial variables X can be obtained via the inverse linear transformation
(17)x=Lz+μ.
Changing variables x→z by using the linear transformation (17) in the integral (10) and taking into account that the Jacobian matrix for this transformation is
(18)∂x∂z=∂(x1,…,xd)∂(z1,…,zd)=L
gives for the integral (10):
(19)E[g]=|det(L)|∫ℝd‍g(Lz+μ)p(Lz+μ)dz.
Taking into account that |det(C)|=|det(LLT)|=|det(L)|2 and substituting the explicit formula for p(x[z]) after the change of variables, one obtains
(20)E[g]=(λ2π)d/2∫ℝd‍g(Lz+μ)exp[-12λzTz]dz.
The Laplace approximation formula for  (11) can be found, for example, in [18]. Applying this formula to integral (20) leads, after a few simplifications, to the asymptotic series:
(21)E[g]~∑m=0∞‍cmλ-m,
where
(22)cm=12m×∑p1+⋯+pd=2m(pkalleven)‍∂2m∂z1p1⋯∂zdpdg(Lz+μ)|z=0×∏k=1d‍1(pk/2)!.
If the expansion is truncated at m=0, the obtained formula is usually referred to as the Laplace approximation with c0=g(μ). For many applications this is sufficient, but not for ours: as will be shown below, the first nontrivial term in the asymptotic expansion for integrals (6)–(9) is linear with respect to (λ-1). Therefore, we retain the first two terms in the asymptotic expansion:
(23)E[g]=g(μ)+12λ∑k=1d‍∂2∂zk2g(Lz+μ)|z=0+O(λ-2).
By introducing the d-dimensional gradient and the Laplace operators, defined as:
(24)∇z=(∂∂z1,…,∂∂zd)T,∇z2=∑k=1d‍∂2∂zk2,
respectively, one can write  (23) in a compact form:
(25)E[g]=g(μ)+12λ∇z2g(Lz+μ)|z=0+O(λ-2).
The subscript next to the nabla operator (e.g., z next to ∇z) will be used to indicate explicitly the variables with respect to which partial derivatives are taken. Asymptotic formulas (23)–(25) can now be applied for the estimation of the integrals (6)–(9) used in the calculation of sensitivity indices.


### 2.3. Calculation of Sensitivity Indices

#### 2.3.1. Variance

In order to calculate parameterD, given by  (7), one can use the well-known statistical formula for the variance:
(26)D=E[f2(x)]-(E[f(x)])2,
where f∅=E[f(x)] is substituted. As one can see from  (26) two integrations are needed: for f(x) and f2(x). We will approximate these integrals independently to the second order, that is, O(λ-2). The application of the approximation formula (25) to the first term in  (26) yields
(27)E[f2(x)]=f2(μ)+12λ∇z2f2(Lz+μ)|z=0+O(λ-2)=f2(μ)+1λ[∇zf(Lz+μ)·∇zf(Lz+μ)+f(μ)∇z2f(Lz+μ)]|z=0+O(λ-2).
Note, that the vector algebra formula
(28)∇2(φψ)=∇·∇(φψ)=2(∇φ)·(∇ψ)+ψ∇2φ+φ∇2ψ
for the Laplacian of the product of two scalar fields, φ and ψ was used in the above derivation.The approximation formula forf∅ can be obtained by a simple substitution of f(x) instead of g(x) into (25), that is,
(29)f∅=E[f(x)]=f(μ)+12λ∇z2f(Lz+μ)|z=0+O(λ-2).
Keeping the terms up to the second order with respect to λ-1 in the asymptotic approximation of f∅2 gives
(30)f∅2=(E[f(x)])2=f2(μ)+f(μ)1λ∇z2f(Lz+μ)|z=0+O(λ-2).
For the sake of the derivation in the next section, let us rewrite  (30) in a different form:
(31)f∅2=(E[f(x)])2=f2(μ)+f(μ)tr(ΣH)+O(λ-2),
where tr(·) denotes the trace of a matrix and H is the Hessian matrix for the function f(x) at a point x=μ with the matrix elements defined as
(32)Hij=∂2f(x)∂xi∂xj|x=μfori,j=1,2,…,d.
After substitution of (27) and (30) into  (26), one obtains the asymptotic approximation to the variance of f:
(33)D=1λ∇zf(Lz+μ)·∇zf(Lz+μ)|z=0+O(λ-2).
Returning back to the initial variables X and taking into account that the gradient operator transforms as
(34)∇zf(Lz+μ)=LT∇xf(x),
where ∇xf(x) is the gradient of f with respect to the “initial” variables Xi (where i=1,2,…,d),  (33) becomes
(35)D=1λ∇xTf(x)LLT∇xf(x)|x=μ+O(λ-2).
Furthermore, taking into account properties (13) and (15), one obtains:
(36)D=∇xTf(x)Σ∇xf(x)|x=μ+O(λ-2).Finally, by introducing the sensitivity vectors=∇xf(x)|x=μ the asymptotic formula for D can be written as:
(37)D=sTΣs+O(λ-2).


#### 2.3.2. Numerator in the Formula for MSI

Formula (8) contains the integrand, which depends on variables from 𝒳u, 𝒳v, and 𝒳v′. Let us introduce the extended vector of input variables:
(38)X~=[XuXvXv′].
Note that the tilde symbol will be used for extended variables. By the assumptions discussed previously Xu, Xv, and Xv′ are statistically independent, hence p(x~)=p(xu)p(xv)p(xv′). Therefore, the random vector X~ has the mean μ~ and the covariance matrix Σ~, which can be written in the following block form:
(39)μ~=[μuμvμv],Σ~=[Σuu000Σvv000Σvv].
Let us introduce auxiliary functions: φ(x~)≡f(xu,xv) and ψ(x~)≡f(xu,xv′) for convenience. Applying formula (28) to g(x~)≡φ(x~)ψ(x~) and following the path used in the previous section, one obtains the asymptotic approximation of the first term in  (8):
(40)E[φ(x~)ψ(x~)]=φ(x~)ψ(x~)|x~=μ~+∇x~Tφ(x~)Σ~∇x~ψ(x~)|x~=μ~+12[φ(x~)tr(Σ~H~)+ψ(x~)tr(Σ~H~′)]|x~=μ~+O(λ-2),
where H~ and H~′ are the Hessian matrices for φ(x~) and ψ(x~), respectively, at x~=μ~. As follows from the definitions of functions φ(x~) and ψ(x~), these matrices have the following block structure:
(41)H~=[HuuHuv0HvuHvv0000],H~′=[Huu0Huv000Hvu0Hvv],
where the blocks are formally defined as
(42)Huu=∂2f(x)∂xu∂xu,Hvv=∂2f(x)∂xv∂xv,Huv=Hvu=∂2f(x)∂xv∂xu.
Similarly, gradients ∇x~φ(x~)|x~=μ~ and ∇x~ψ(x~)|x~=μ~ have the following matrix structure:
(43)∇x~φ(x~)|x~=μ~=[susv0],∇x~ψ(x~)|x~=μ~=[su0sv],
where su=∇xuf(x)|x=μ and sv=∇xvf(x)|x=μ are sensitivity vectors with respect to Xu and Xv. Taking into account (41) and (43), the trace operators in  (40) can be modified in the following way: tr(Σ~H~)=tr(Σ~H~′)=tr(ΣuuHuu)+tr(ΣvvHvv)=tr(ΣH). From the definitions of φ(x~) and ψ(x~) it follows that φ(μ~)=ψ(μ~)=f(μ), hence  (43) becomes
(44)E[φ(x~)ψ(x~)]=f2(μ)+suTΣuusu+f(μ)tr(ΣH)+O(λ-2).
After subtracting  (31) from  (44) one obtains the following approximation for the factor DXu in the formula for the Main Sensitivity Index:
(45)DXu=suTΣuusu+O(λ-2).


#### 2.3.3. Numerator in the Formula for TSI

The calculation ofDXutot requires the calculation of integral (9). In this section let us introduce the extended variable vector as the column vector
(46)X~=[XuXvXu′].
Vector X~ has the mean and covariance matrix that have the following block representation:
(47)μ~=[μuμvμu],Σ~=[Σuu000Σvv000Σuu].
Introducing, for convenience, the auxiliary functions: φ(x~)≡f(xu,xv), χ(x~)≡f(xu′,xv) and ω(x~)≡φ(x~)-χ(x~), one can write function g(x~) in  (10) as
(48)g(x~)=12ω2(x~).
The Laplace approximation for DXutot=E[g(x~)] can be obtained by substituting function (48) into formula (25), which gives:
(49)DXutot=12ω2(x~)|x~=μ~+14λ~∇z~2ω2(L~z~+μ~)|z~=0+O(λ-2).
From the definition of functions φ(x~) and χ(x~) it follows that ω(μ~)=0 and the first term in  (49) vanishes. The second term can be transformed in the following way:
(50)DXutot=14λ~∇z~2ω2(L~z~+μ~)|z~=0+O(λ-2)=12λ~[∇z~ω(L~z~+μ~)·∇z~ω(L~z~+μ~)+ω(μ~)·∇z~2ω(L~z~+μ~)]|z~=0+O(λ-2)=12∇x~Tω(x~)Σ~∇x~ω(x~)|x~=μ~+O(λ-2).
Taking into account the structure of the extended covariance matrix, given by  (47), and the structure of the gradients in  (50), given by
(51)∇x~ω(x~)|x~=μ~=[-su0-su],
one obtains the following formula:
(52)DXutot=suTΣuusu+O(λ-2).


#### 2.3.4. Asymptotic Formulas for the Global Sensitivity Indices

In the previous sections we have obtained approximation formulas for the integrals involved in the calculation of the global sensitivity indices. As we have demonstrated, these integrals can be approximated with the asymptotic series∑m=1∞‍cmλ-m, where the coefficients cm depend on the derivatives of the model function with respect to input variables. Based on this result, the asymptotic formula for the sensitivity indices, SXu and SXutot can be written as a ratio:
(53)∑n=1∞anλ-n∑k=1∞bkλ-k=a1b1+(a2b1-a1b2)b12λ-1+O(λ-2).
The first two terms of the Maclaurin series with respect to the small parameter λ-1 are written explicitly on the right-hand side of  (53). Hence the first nontrivial term in the approximation of the sensitivity indices is the ratio of the first nonvanishing term in the integral approximations obtained above. Thus, the approximation for MSI can be obtained by dividing  (45) by  (37):
(54)SXu=suTΣuususTΣs+O(λ-1).
Similarly, the approximation for the TSI is
(55)SXutot=suTΣuususTΣs+O(λ-1).


## 2.1. Variance-Based Global Sensitivity Analysis

Let𝒳={X1,X2,…,Xd} be a set of continuous random variables, 𝒳u be its arbitrary subset and 𝒳v be a subset of complimentary variables, that is, 𝒳=𝒳u∪𝒳v. The random variables can be gathered into the corresponding column vectors: Xu, Xv, and X=(Xu,Xv). Following the arguments given in [5, 15, 16], we assume that random variables from 𝒳u and 𝒳v are mutually statistically independent, hence their joint probability function is the product of the marginal probability density functions:
(2)p(x)=p(xu)p(xv).
Note that we will use, as is the rule in statistics, a capital letter (e.g., Xi) to denote a random variable and the corresponding lowercase letter (xi) to denote its value (realizations). Let us assume that the random variables Xi(i=1,2,…,d) are distributed according to the normal law with known means and covariances. The multivariate normal distribution is characterized by the probability density function [17]
(3)p(x)=1(2π)d/2det(Σ)1/2exp[-12(x-μ)TΣ-1(x-μ)],
where X is the column vector of the input random variables, μ=E[X] is the column vector of their expected values, Σ=E[(x-μ)(x-μ)T] is the covariance matrix and E denotes the mathematical expectation and the symbol “T” denotes the operation of transposing a row to a column.Consider a model(4)Y=f(X1,X2,…,Xd)=f(X),
where the function f:ℝd→ℝ is generally nonlinear and may result from numerical calculations. Xi are called inputs and Y is called the output or response. The response Y=f(X) is a random variable itself, with the expected value f∅=E[f(X)] and the variance D=E[(f(X)-f∅)2]. The variance of the output, D, can be used to characterize the uncertainty in the response due to the uncertainty in the input, or in other words, the uncertainty propagation. Note that D can be calculated, regardless of whether the normally distributed inputs are correlated or not. Sensitivity indices characterize the contribution of a subset 𝒳u (and, correspondingly, vector Xu) to the response variance and they are meaningful only if subset 𝒳u and its complementary subset 𝒳v are statistically independent. Two types of global sensitivity indices, with respect to a subset 𝒳u, can be introduced:(1)
SXu, which is called the main sensitivity index (MSI) and it represents the effect due to Xu only;(2)
SXutot, which is called the total sensitivity index (TSI) and it represents the contribution to the variance of Xu along with all the interactions of this variable with other variables.The global sensitivity indices can be calculated using the ratios [1, 8]
(5)SXutot=DXutotD,SXu=DXuD,
where parameters D, DXu, and DXutot can be calculated through the following multidimensional integrals [12, 15, 16]:
(6)f∅=∫ℝd‍f(x)p(x)dx,(7)D=∫ℝd‍f2(x)p(x)dx-f∅2,(8)DXu=∫ℝ2d-dim(Xu)‍f(x)f(xu,xv′)p(x)p(xv′)dxdxv′-f∅2,(9)DXutot=12∫ℝd+dim(Xu)‍[f(x)-f(xu′,xv)]2p(x)p(xu′)dxdxu′,
which are a generalization of the classical formulas given by Sobol’ [1, 8]. Here the prime symbol over a variable (e.g., as in xu′) means that this variable has to be sampled independently from the corresponding marginal distribution (p(xu′) in this case) of its unprimed analogue and dim(·) is the number of elements in a vector, that is, the dimensionality of the vector. Calculation of the sensitivity indices requires the calculation of integrals (6)–(9), which can be written in the following general form:
(10)E[g]=∫ℝdeff‍g(x~)p(x~)dx~,
where the integration operator E[·] denotes the mathematical expectation with respect to extended random variables X~, g(x~) represents a function being integrated and deff=dim(X~) is the effective dimensionality of the integral. For instance, in integral (8) function g(x~) represents [f(x)f(xu,xv′)], x~ represents (xu,xv,xv′), and the effective dimensionality is deff=2d-dim(Xu). Similarly, in the integral (6) function g(x~) represents f(x), x~=x and deff=d.


## 2.2. Laplace Approximation of Integrals

In the limit when elements of the covariance matrix are small in some norm, one can use for the calculation of integrals (6)–(9) the asymptotic analysis technique known as the Laplace approximation. The Laplace approximation is an approximation to integrals in the following form:
(11)I(λ)=∫Ω‍h(x)exp[λR(x)]dx,
as λ→∞ with λ>0. Here Ω is a measurable subset of ℝd and R(x) is a real-valued function of x. When applying the Laplace approximation formula, it is assumed that R(x) is an infinity differentiable function with its maximum value at a single point x0 in the interior of Ω (see [18] for details) and h(x) is a sufficiently smooth function of its argument.Integrals in the form (10), used for the calculation of sensitivity indices, fulfill the above requirements for the applicability of the Laplace approximation formula if one considers x0=μ and assumes that elements of the covariance matrix Σ are small. Let us assume that the covariance matrix Σ is small in some norm, that is, in the limit ∥Σ∥→0. For the sake of the asymptotic analysis, let us introduce the big parameter λ as the inverse absolute value of the determinant of the covariance matrix
(12)λ=|det(Σ)|-1.
We will assume that the covariance matrix is strictly positive definite, hence the absolute value sign in (12) may be omitted. Factoring out the coefficient λ allows us to introduce the normalized covariance matrix
(13)C=Σ|det(Σ)|-1=λΣ.
Taking into account that det(C)=det(λΣ)=λddet(Σ), the probability density function (3) can be written in a form suitable for the Laplace approximation:
(14)p(x)=1|det(C)|1/2(λ2π)d/2exp[-λ12(x-μ)TC-1(x-μ)].
The covariance matrix Σ is real-valued, symmetric, and positively definite, hence matrix C allows the Cholesky factorization:
(15)C=LLT,
where L is a lower triangular matrix with strictly positive diagonal entries. The exponent in the joint probability density function can then be transformed in the following way:
(16)λ12(x-μ)TC-1(x-μ)=λ12zTz,
where new random variables Zi(i=1,…,d) are introduced by the linear transformation z=L-1(x-μ). The initial variables X can be obtained via the inverse linear transformation
(17)x=Lz+μ.
Changing variables x→z by using the linear transformation (17) in the integral (10) and taking into account that the Jacobian matrix for this transformation is
(18)∂x∂z=∂(x1,…,xd)∂(z1,…,zd)=L
gives for the integral (10):
(19)E[g]=|det(L)|∫ℝd‍g(Lz+μ)p(Lz+μ)dz.
Taking into account that |det(C)|=|det(LLT)|=|det(L)|2 and substituting the explicit formula for p(x[z]) after the change of variables, one obtains
(20)E[g]=(λ2π)d/2∫ℝd‍g(Lz+μ)exp[-12λzTz]dz.
The Laplace approximation formula for  (11) can be found, for example, in [18]. Applying this formula to integral (20) leads, after a few simplifications, to the asymptotic series:
(21)E[g]~∑m=0∞‍cmλ-m,
where
(22)cm=12m×∑p1+⋯+pd=2m(pkalleven)‍∂2m∂z1p1⋯∂zdpdg(Lz+μ)|z=0×∏k=1d‍1(pk/2)!.
If the expansion is truncated at m=0, the obtained formula is usually referred to as the Laplace approximation with c0=g(μ). For many applications this is sufficient, but not for ours: as will be shown below, the first nontrivial term in the asymptotic expansion for integrals (6)–(9) is linear with respect to (λ-1). Therefore, we retain the first two terms in the asymptotic expansion:
(23)E[g]=g(μ)+12λ∑k=1d‍∂2∂zk2g(Lz+μ)|z=0+O(λ-2).
By introducing the d-dimensional gradient and the Laplace operators, defined as:
(24)∇z=(∂∂z1,…,∂∂zd)T,∇z2=∑k=1d‍∂2∂zk2,
respectively, one can write  (23) in a compact form:
(25)E[g]=g(μ)+12λ∇z2g(Lz+μ)|z=0+O(λ-2).
The subscript next to the nabla operator (e.g., z next to ∇z) will be used to indicate explicitly the variables with respect to which partial derivatives are taken. Asymptotic formulas (23)–(25) can now be applied for the estimation of the integrals (6)–(9) used in the calculation of sensitivity indices.


## 2.3. Calculation of Sensitivity Indices

### 2.3.1. Variance

In order to calculate parameterD, given by  (7), one can use the well-known statistical formula for the variance:
(26)D=E[f2(x)]-(E[f(x)])2,
where f∅=E[f(x)] is substituted. As one can see from  (26) two integrations are needed: for f(x) and f2(x). We will approximate these integrals independently to the second order, that is, O(λ-2). The application of the approximation formula (25) to the first term in  (26) yields
(27)E[f2(x)]=f2(μ)+12λ∇z2f2(Lz+μ)|z=0+O(λ-2)=f2(μ)+1λ[∇zf(Lz+μ)·∇zf(Lz+μ)+f(μ)∇z2f(Lz+μ)]|z=0+O(λ-2).
Note, that the vector algebra formula
(28)∇2(φψ)=∇·∇(φψ)=2(∇φ)·(∇ψ)+ψ∇2φ+φ∇2ψ
for the Laplacian of the product of two scalar fields, φ and ψ was used in the above derivation.The approximation formula forf∅ can be obtained by a simple substitution of f(x) instead of g(x) into (25), that is,
(29)f∅=E[f(x)]=f(μ)+12λ∇z2f(Lz+μ)|z=0+O(λ-2).
Keeping the terms up to the second order with respect to λ-1 in the asymptotic approximation of f∅2 gives
(30)f∅2=(E[f(x)])2=f2(μ)+f(μ)1λ∇z2f(Lz+μ)|z=0+O(λ-2).
For the sake of the derivation in the next section, let us rewrite  (30) in a different form:
(31)f∅2=(E[f(x)])2=f2(μ)+f(μ)tr(ΣH)+O(λ-2),
where tr(·) denotes the trace of a matrix and H is the Hessian matrix for the function f(x) at a point x=μ with the matrix elements defined as
(32)Hij=∂2f(x)∂xi∂xj|x=μfori,j=1,2,…,d.
After substitution of (27) and (30) into  (26), one obtains the asymptotic approximation to the variance of f:
(33)D=1λ∇zf(Lz+μ)·∇zf(Lz+μ)|z=0+O(λ-2).
Returning back to the initial variables X and taking into account that the gradient operator transforms as
(34)∇zf(Lz+μ)=LT∇xf(x),
where ∇xf(x) is the gradient of f with respect to the “initial” variables Xi (where i=1,2,…,d),  (33) becomes
(35)D=1λ∇xTf(x)LLT∇xf(x)|x=μ+O(λ-2).
Furthermore, taking into account properties (13) and (15), one obtains:
(36)D=∇xTf(x)Σ∇xf(x)|x=μ+O(λ-2).Finally, by introducing the sensitivity vectors=∇xf(x)|x=μ the asymptotic formula for D can be written as:
(37)D=sTΣs+O(λ-2).


### 2.3.2. Numerator in the Formula for MSI

Formula (8) contains the integrand, which depends on variables from 𝒳u, 𝒳v, and 𝒳v′. Let us introduce the extended vector of input variables:
(38)X~=[XuXvXv′].
Note that the tilde symbol will be used for extended variables. By the assumptions discussed previously Xu, Xv, and Xv′ are statistically independent, hence p(x~)=p(xu)p(xv)p(xv′). Therefore, the random vector X~ has the mean μ~ and the covariance matrix Σ~, which can be written in the following block form:
(39)μ~=[μuμvμv],Σ~=[Σuu000Σvv000Σvv].
Let us introduce auxiliary functions: φ(x~)≡f(xu,xv) and ψ(x~)≡f(xu,xv′) for convenience. Applying formula (28) to g(x~)≡φ(x~)ψ(x~) and following the path used in the previous section, one obtains the asymptotic approximation of the first term in  (8):
(40)E[φ(x~)ψ(x~)]=φ(x~)ψ(x~)|x~=μ~+∇x~Tφ(x~)Σ~∇x~ψ(x~)|x~=μ~+12[φ(x~)tr(Σ~H~)+ψ(x~)tr(Σ~H~′)]|x~=μ~+O(λ-2),
where H~ and H~′ are the Hessian matrices for φ(x~) and ψ(x~), respectively, at x~=μ~. As follows from the definitions of functions φ(x~) and ψ(x~), these matrices have the following block structure:
(41)H~=[HuuHuv0HvuHvv0000],H~′=[Huu0Huv000Hvu0Hvv],
where the blocks are formally defined as
(42)Huu=∂2f(x)∂xu∂xu,Hvv=∂2f(x)∂xv∂xv,Huv=Hvu=∂2f(x)∂xv∂xu.
Similarly, gradients ∇x~φ(x~)|x~=μ~ and ∇x~ψ(x~)|x~=μ~ have the following matrix structure:
(43)∇x~φ(x~)|x~=μ~=[susv0],∇x~ψ(x~)|x~=μ~=[su0sv],
where su=∇xuf(x)|x=μ and sv=∇xvf(x)|x=μ are sensitivity vectors with respect to Xu and Xv. Taking into account (41) and (43), the trace operators in  (40) can be modified in the following way: tr(Σ~H~)=tr(Σ~H~′)=tr(ΣuuHuu)+tr(ΣvvHvv)=tr(ΣH). From the definitions of φ(x~) and ψ(x~) it follows that φ(μ~)=ψ(μ~)=f(μ), hence  (43) becomes
(44)E[φ(x~)ψ(x~)]=f2(μ)+suTΣuusu+f(μ)tr(ΣH)+O(λ-2).
After subtracting  (31) from  (44) one obtains the following approximation for the factor DXu in the formula for the Main Sensitivity Index:
(45)DXu=suTΣuusu+O(λ-2).


### 2.3.3. Numerator in the Formula for TSI

The calculation ofDXutot requires the calculation of integral (9). In this section let us introduce the extended variable vector as the column vector
(46)X~=[XuXvXu′].
Vector X~ has the mean and covariance matrix that have the following block representation:
(47)μ~=[μuμvμu],Σ~=[Σuu000Σvv000Σuu].
Introducing, for convenience, the auxiliary functions: φ(x~)≡f(xu,xv), χ(x~)≡f(xu′,xv) and ω(x~)≡φ(x~)-χ(x~), one can write function g(x~) in  (10) as
(48)g(x~)=12ω2(x~).
The Laplace approximation for DXutot=E[g(x~)] can be obtained by substituting function (48) into formula (25), which gives:
(49)DXutot=12ω2(x~)|x~=μ~+14λ~∇z~2ω2(L~z~+μ~)|z~=0+O(λ-2).
From the definition of functions φ(x~) and χ(x~) it follows that ω(μ~)=0 and the first term in  (49) vanishes. The second term can be transformed in the following way:
(50)DXutot=14λ~∇z~2ω2(L~z~+μ~)|z~=0+O(λ-2)=12λ~[∇z~ω(L~z~+μ~)·∇z~ω(L~z~+μ~)+ω(μ~)·∇z~2ω(L~z~+μ~)]|z~=0+O(λ-2)=12∇x~Tω(x~)Σ~∇x~ω(x~)|x~=μ~+O(λ-2).
Taking into account the structure of the extended covariance matrix, given by  (47), and the structure of the gradients in  (50), given by
(51)∇x~ω(x~)|x~=μ~=[-su0-su],
one obtains the following formula:
(52)DXutot=suTΣuusu+O(λ-2).


### 2.3.4. Asymptotic Formulas for the Global Sensitivity Indices

In the previous sections we have obtained approximation formulas for the integrals involved in the calculation of the global sensitivity indices. As we have demonstrated, these integrals can be approximated with the asymptotic series∑m=1∞‍cmλ-m, where the coefficients cm depend on the derivatives of the model function with respect to input variables. Based on this result, the asymptotic formula for the sensitivity indices, SXu and SXutot can be written as a ratio:
(53)∑n=1∞anλ-n∑k=1∞bkλ-k=a1b1+(a2b1-a1b2)b12λ-1+O(λ-2).
The first two terms of the Maclaurin series with respect to the small parameter λ-1 are written explicitly on the right-hand side of  (53). Hence the first nontrivial term in the approximation of the sensitivity indices is the ratio of the first nonvanishing term in the integral approximations obtained above. Thus, the approximation for MSI can be obtained by dividing  (45) by  (37):
(54)SXu=suTΣuususTΣs+O(λ-1).
Similarly, the approximation for the TSI is
(55)SXutot=suTΣuususTΣs+O(λ-1).


## 2.3.1. Variance

In order to calculate parameterD, given by  (7), one can use the well-known statistical formula for the variance:
(26)D=E[f2(x)]-(E[f(x)])2,
where f∅=E[f(x)] is substituted. As one can see from  (26) two integrations are needed: for f(x) and f2(x). We will approximate these integrals independently to the second order, that is, O(λ-2). The application of the approximation formula (25) to the first term in  (26) yields
(27)E[f2(x)]=f2(μ)+12λ∇z2f2(Lz+μ)|z=0+O(λ-2)=f2(μ)+1λ[∇zf(Lz+μ)·∇zf(Lz+μ)+f(μ)∇z2f(Lz+μ)]|z=0+O(λ-2).
Note, that the vector algebra formula
(28)∇2(φψ)=∇·∇(φψ)=2(∇φ)·(∇ψ)+ψ∇2φ+φ∇2ψ
for the Laplacian of the product of two scalar fields, φ and ψ was used in the above derivation.The approximation formula forf∅ can be obtained by a simple substitution of f(x) instead of g(x) into (25), that is,
(29)f∅=E[f(x)]=f(μ)+12λ∇z2f(Lz+μ)|z=0+O(λ-2).
Keeping the terms up to the second order with respect to λ-1 in the asymptotic approximation of f∅2 gives
(30)f∅2=(E[f(x)])2=f2(μ)+f(μ)1λ∇z2f(Lz+μ)|z=0+O(λ-2).
For the sake of the derivation in the next section, let us rewrite  (30) in a different form:
(31)f∅2=(E[f(x)])2=f2(μ)+f(μ)tr(ΣH)+O(λ-2),
where tr(·) denotes the trace of a matrix and H is the Hessian matrix for the function f(x) at a point x=μ with the matrix elements defined as
(32)Hij=∂2f(x)∂xi∂xj|x=μfori,j=1,2,…,d.
After substitution of (27) and (30) into  (26), one obtains the asymptotic approximation to the variance of f:
(33)D=1λ∇zf(Lz+μ)·∇zf(Lz+μ)|z=0+O(λ-2).
Returning back to the initial variables X and taking into account that the gradient operator transforms as
(34)∇zf(Lz+μ)=LT∇xf(x),
where ∇xf(x) is the gradient of f with respect to the “initial” variables Xi (where i=1,2,…,d),  (33) becomes
(35)D=1λ∇xTf(x)LLT∇xf(x)|x=μ+O(λ-2).
Furthermore, taking into account properties (13) and (15), one obtains:
(36)D=∇xTf(x)Σ∇xf(x)|x=μ+O(λ-2).Finally, by introducing the sensitivity vectors=∇xf(x)|x=μ the asymptotic formula for D can be written as:
(37)D=sTΣs+O(λ-2).


## 2.3.2. Numerator in the Formula for MSI

Formula (8) contains the integrand, which depends on variables from 𝒳u, 𝒳v, and 𝒳v′. Let us introduce the extended vector of input variables:
(38)X~=[XuXvXv′].
Note that the tilde symbol will be used for extended variables. By the assumptions discussed previously Xu, Xv, and Xv′ are statistically independent, hence p(x~)=p(xu)p(xv)p(xv′). Therefore, the random vector X~ has the mean μ~ and the covariance matrix Σ~, which can be written in the following block form:
(39)μ~=[μuμvμv],Σ~=[Σuu000Σvv000Σvv].
Let us introduce auxiliary functions: φ(x~)≡f(xu,xv) and ψ(x~)≡f(xu,xv′) for convenience. Applying formula (28) to g(x~)≡φ(x~)ψ(x~) and following the path used in the previous section, one obtains the asymptotic approximation of the first term in  (8):
(40)E[φ(x~)ψ(x~)]=φ(x~)ψ(x~)|x~=μ~+∇x~Tφ(x~)Σ~∇x~ψ(x~)|x~=μ~+12[φ(x~)tr(Σ~H~)+ψ(x~)tr(Σ~H~′)]|x~=μ~+O(λ-2),
where H~ and H~′ are the Hessian matrices for φ(x~) and ψ(x~), respectively, at x~=μ~. As follows from the definitions of functions φ(x~) and ψ(x~), these matrices have the following block structure:
(41)H~=[HuuHuv0HvuHvv0000],H~′=[Huu0Huv000Hvu0Hvv],
where the blocks are formally defined as
(42)Huu=∂2f(x)∂xu∂xu,Hvv=∂2f(x)∂xv∂xv,Huv=Hvu=∂2f(x)∂xv∂xu.
Similarly, gradients ∇x~φ(x~)|x~=μ~ and ∇x~ψ(x~)|x~=μ~ have the following matrix structure:
(43)∇x~φ(x~)|x~=μ~=[susv0],∇x~ψ(x~)|x~=μ~=[su0sv],
where su=∇xuf(x)|x=μ and sv=∇xvf(x)|x=μ are sensitivity vectors with respect to Xu and Xv. Taking into account (41) and (43), the trace operators in  (40) can be modified in the following way: tr(Σ~H~)=tr(Σ~H~′)=tr(ΣuuHuu)+tr(ΣvvHvv)=tr(ΣH). From the definitions of φ(x~) and ψ(x~) it follows that φ(μ~)=ψ(μ~)=f(μ), hence  (43) becomes
(44)E[φ(x~)ψ(x~)]=f2(μ)+suTΣuusu+f(μ)tr(ΣH)+O(λ-2).
After subtracting  (31) from  (44) one obtains the following approximation for the factor DXu in the formula for the Main Sensitivity Index:
(45)DXu=suTΣuusu+O(λ-2).


## 2.3.3. Numerator in the Formula for TSI

The calculation ofDXutot requires the calculation of integral (9). In this section let us introduce the extended variable vector as the column vector
(46)X~=[XuXvXu′].
Vector X~ has the mean and covariance matrix that have the following block representation:
(47)μ~=[μuμvμu],Σ~=[Σuu000Σvv000Σuu].
Introducing, for convenience, the auxiliary functions: φ(x~)≡f(xu,xv), χ(x~)≡f(xu′,xv) and ω(x~)≡φ(x~)-χ(x~), one can write function g(x~) in  (10) as
(48)g(x~)=12ω2(x~).
The Laplace approximation for DXutot=E[g(x~)] can be obtained by substituting function (48) into formula (25), which gives:
(49)DXutot=12ω2(x~)|x~=μ~+14λ~∇z~2ω2(L~z~+μ~)|z~=0+O(λ-2).
From the definition of functions φ(x~) and χ(x~) it follows that ω(μ~)=0 and the first term in  (49) vanishes. The second term can be transformed in the following way:
(50)DXutot=14λ~∇z~2ω2(L~z~+μ~)|z~=0+O(λ-2)=12λ~[∇z~ω(L~z~+μ~)·∇z~ω(L~z~+μ~)+ω(μ~)·∇z~2ω(L~z~+μ~)]|z~=0+O(λ-2)=12∇x~Tω(x~)Σ~∇x~ω(x~)|x~=μ~+O(λ-2).
Taking into account the structure of the extended covariance matrix, given by  (47), and the structure of the gradients in  (50), given by
(51)∇x~ω(x~)|x~=μ~=[-su0-su],
one obtains the following formula:
(52)DXutot=suTΣuusu+O(λ-2).


## 2.3.4. Asymptotic Formulas for the Global Sensitivity Indices

In the previous sections we have obtained approximation formulas for the integrals involved in the calculation of the global sensitivity indices. As we have demonstrated, these integrals can be approximated with the asymptotic series∑m=1∞‍cmλ-m, where the coefficients cm depend on the derivatives of the model function with respect to input variables. Based on this result, the asymptotic formula for the sensitivity indices, SXu and SXutot can be written as a ratio:
(53)∑n=1∞anλ-n∑k=1∞bkλ-k=a1b1+(a2b1-a1b2)b12λ-1+O(λ-2).
The first two terms of the Maclaurin series with respect to the small parameter λ-1 are written explicitly on the right-hand side of  (53). Hence the first nontrivial term in the approximation of the sensitivity indices is the ratio of the first nonvanishing term in the integral approximations obtained above. Thus, the approximation for MSI can be obtained by dividing  (45) by  (37):
(54)SXu=suTΣuususTΣs+O(λ-1).
Similarly, the approximation for the TSI is
(55)SXutot=suTΣuususTΣs+O(λ-1).


## 3. Discussion

Formula (37) relates the input uncertainties, given by the covariance matrix Σ, to the output uncertainty, characterized by the variance D and, hence, describes the propagation of error. Contrary to what one would expect from  (23), the approximation formula (36) involves first-order partial derivatives of f(x) and not zeroth and second-order derivatives.As one can see, the structure of the first nonvanishing term in this expression corresponds to the famous sandwich formula for error propagation used in the derivative-based LSA [3, 4]. Note, that in this work sensitivity vector components are defined as si=(∂f(x)/∂xi)x=μ, whereas in literature dimensionless sensitivity coefficients si=(μi/f(μ))(∂f(x)/∂xi)x=μ are often used.Formulas (45) and (52) have the same sandwich structure as  (37), but contains only the sensitivities and uncertainties with respect to the input variables Xu. Moreover, factor DXu matches factor DXutot under conditions used in the approximation and the corresponding sensitivity indices coincide. As a result, sensitivity indices given by (54) and (55) do not differ in the first nonvanishing term in the approximation. The difference between SXu and SXutot may be of the order O(λ-1). Formulas (54) and (55) involve, besides elements of the covariance matrix, only gradients of the model function, which are proportional to the coefficients of the linear expansion of the model at x=μ. This indicates that, in the considered limit, the model is approximately linear. As a consequence, interactions are negligible and SXu coincides with SXu that is reflected by results (54) and (55).The remaining question of this analysis is the criterion for the applicability of the asymptotic results. We could not derive this criterion directly for the sensitivity indices. The indirect indication can obtained by using  (27):
(56)f2(μ)≫sTΣs+f(μ)tr(ΣH),
by taking into account that the calculation of the sensitivity index involves the calculation of integrals for f2(x).


## 4. Example

Adetula and Bokov [15, 16] discussed a test problem, where the model f:ℝ7→ℝ represents the two-group infinite multiplication factor described by the analytical formula [19]:
(57)f(X1,…,X7)=X5X1+X3+X7+X6X7(X2+X4)(X1+X3+X7),
where variables Xi with i=1,2,…,7 represent the following macroscopic cross-sections: X1 is the fast capture, X2 is the thermal capture, X3 is the fast fission, X4 is the thermal fission, X5 is the fast neutron production, X6 is the thermal neutron production, and X7 is the fast removal. The mean values of the cross-sections and covariance matrix used were based on the data given in [20]:
(58)μ=10-2×[0.53362.6930.191242.84380.49206.9292.063]T,Σ=10-8×[0.420.22-0.05-0.020000.222.14-0.02-0.65000-0.05-0.020.020.01000-0.02-0.650.010.8400000000.230.18000000.189.6400000005.28],
where the cross-sections are given in cm-1 and the variances are given in cm-2. The test covariance matrix has three uncorrelated blocks, corresponding to subsets: {X1,X2,X3,X4}, {X4,X5}, and {X7}, for which the calculation of sensitivity indices is not meaningless.Values of the global sensitivity indices, which are calculated using the asymptotic formulas obtained in our work and the values as reported in [15, 16] are presented in Table 1. The output uncertainty, characterized by the variance (D) or the standard deviation (D) calculated with these two methods are also given in Table 1. The values taken from [15, 16] were calculated with the sparse grid quadrature, as it has been found to be the most accurate amongst all the applied quadratures. (Other quadratures were Monte Carlo and quasi-Monte Carlo. Values obtained with all three quadratures were in a good agreement, this indicates that the quadratures have converged and the result obtained may be taken as the reference.) As one can see from Table 1, the results for sensitivity indices and uncertainties are in good agreement and differ for only one case in the fourth significant digit for sensitivity indices and in the third significant digit for the variance. It is worth mentioning that this discrepancy is still within the error of numerical integration as reported in [15, 16]. Results given in Table 1 confirm the validity of the asymptotic approximation for the considered test problem. Similar results have been obtained for the second case, as discussed in [15, 16], namely the case of the diagonal covariance matrix. These results also confirm the correctness of the hypothesis given in [15, 16] and discussed in earlier in Section 1 of this paper.Table 1
Sensitivity indices and output uncertainty (both dimensionless) as calculated with the asymptotic formula and multidimensional quadrature.
Parameter
Asymptotic approximation
Sparse grid integration [15, 16]
Sensitivity index for{X1,X2,X3,X4}
0.3453
0.3453
Sensitivity index for{X5,X6}
0.6126
0.6125
Sensitivity index for{X7}
0.0421
0.0421
VarianceD
3.575
·
1
0
-
5
3.567
·
1
0
-
5
Standard deviationD
598 pcm
597 pcm


## 5. Conclusions

The problem of calculating the uncertainty and sensitivity parameters in the framework of the variance-based Global Sensitivity Analysis was addressed. It was assumed that the input variables are normally distributed and block-wise correlated. Under the additional assumption that the covariance matrix is small in some norm, the Laplace approximation technique was applied and the asymptotic approximation of the output variance and the global sensitivity indices was obtained.Our results demonstrate that the first nontrivial term in the asymptotic expansion of the output uncertainty (variance) has the so-called sandwich structure well known from the local sensitivity analysis, that is, it depends linearly on the covariance matrix and quadratically on the partial derivatives of the outputs with respect to the inputs. The dependence of uncertainty on the partial derivatives means that, under the above assumptions, the model can be considered to be approximately linear, and, as a consequence, the interaction (but not correlation) between inputs is negligible. The above conclusion is supported by the asymptotic approximation obtained for sensitivity indices: the expression for the main sensitivity indices coincides with the expression for the total sensitivity indices, thus indicating that interactions between input variables can be neglected.A test problem with realistic values of multigroup cross-sections and their covariances was considered. The values of the global sensitivity indices calculated via asymptotic formula were found to be in good agreement with the values calculated from the exact formulas using numerical integration. This demonstrate that the results obtained may not only be of academic interest but of practical interest as well. The theoretical and numerical results obtained in this work indicate that, from the point of view of numerical calculations, the sparse-grid integration-based method for the uncertainty propagation can be an alternative to the perturbation/derivative methods used in the local sensitivity analysis. Moreover, the results for the uncertainty obtained by these two methods are interchangeable in the limit of small cross-section errors.



---

*Source: 253045-2012-09-19.xml*