Efficient computation of global sensitivity indices using sparse polynomial chaos expansions

https://doi.org/10.1016/j.ress.2010.06.015Get rights and content

Abstract

Global sensitivity analysis aims at quantifying the relative importance of uncertain input variables onto the response of a mathematical model of a physical system. ANOVA-based indices such as the Sobol’ indices are well-known in this context. These indices are usually computed by direct Monte Carlo or quasi-Monte Carlo simulation, which may reveal hardly applicable for computationally demanding industrial models. In the present paper, sparse polynomial chaos (PC) expansions are introduced in order to compute sensitivity indices. An adaptive algorithm allows the analyst to build up a PC-based metamodel that only contains the significant terms whereas the PC coefficients are computed by least-square regression using a computer experimental design. The accuracy of the metamodel is assessed by leave-one-out cross validation. Due to the genuine orthogonality properties of the PC basis, ANOVA-based sensitivity indices are post-processed analytically. This paper also develops a bootstrap technique which eventually yields confidence intervals on the results. The approach is illustrated on various application examples up to 21 stochastic dimensions. Accurate results are obtained at a computational cost 2–3 orders of magnitude smaller than that associated with Monte Carlo simulation.

Introduction

The mathematical model of a physical system can be regarded as a function of a set of input parameters. In many cases, such a mapping does not have an explicit analytical expression due to the complexity of the underlying phenomena. Evaluating the model response often requires to solve a system of differential equations, which can be achieved by employing discretization schemes (e.g. finite difference and finite element methods) that are usually implemented in computer codes. The increasing availability of computational power and complexity of numerical solvers have allowed one to reduce considerably the numerical errors, and to handle more and more realistic behaviours. However, numerical simulation can provide realistic predictions only if the model input parameters are accurately selected. This is seldom the case in practice since the latter are often affected by uncertainty, which may either result from a lack of knowledge (epistemic uncertainty) or be associated to an intrinsic variability (aleatory uncertainty). In this context, quantifying the contribution of each input parameter to the output variability is of major interest. This may not be an easy task in practice though since the relationship between the input and the output variables of a complex model is not straightforward.

Sensitivity analysis allows one to address this problem. Various methods have been investigated in the literature [1], [2], which are typically separated into two groups:

  • local sensitivity analysis, which studies how little variations of the input parameters in the vicinity of given values influence the model response;

  • global sensitivity analysis, which is related with quantifying the output uncertainty due to changes of the input parameters (which are taken singly or in combination with others) over their entire domain of variation.

Many papers have been devoted to the latter topic in the last two decades. A state of the art review is available in [1], which gathers the related methods into two categories:
  • regression-based methods, which exploit the results of the linear regression of the model response onto the input vector. These approaches are useful to measure the effects of the input variables if the model is linear or almost linear. However, they fail to produce satisfactory sensitivity measures in case of significant nonlinearity [3].

  • variance-based methods, which rely upon the decomposition of the response variance as a sum of contributions of each input variables, or combinations thereof. They are known as ANOVA (for ANalysis Of VAriance) techniques in statistics [4]. The Fourier amplitude sensitivity test (FAST) indices [5], [6] and the Sobol’ indices [3], [7], [8], [9] enter this category, see also the review in [10].

Other global sensitivity analysis techniques are available, such as Morris method [11], sampling methods [12] and methods based on experimental designs, e.g. fractional factorial designs [13] and Plackett–Burman designs [14].

Of interest are the Sobol’ sensitivity indices in this paper. These quantities are often computed using crude Monte Carlo simulation, which may make them hardly applicable for computationally demanding models. To overcome this difficulty, one may substitute the model under consideration for an analytical approximation, named metamodel, whose evaluations are inexpensive. The use of various metamodels based on non-parametric regression procedures was investigated in [15], [16], [17]. It is also worth mentioning specific developments of the so-called Gaussian process modelling approach [18], [19]. As an alternative, the metamodelling of the functions arising from the ANOVA decomposition rather than the model function was proposed in [20]. The approach belongs to the family of the high dimensional model representation (HDMR) methods, which consist in neglecting the ANOVA functionals associated with high-order interactions. The use of orthogonal polynomial bases has received much interest in this context [21], [22], [23]. This choice is relevant since the Sobol’ indices can be then computed exactly from algebraic operations on the coefficients of the polynomial expansions.

The present paper is also focused on representations onto polynomial bases, but in a more classical metamodelling framework since one approximates directly the model function. The approach devised herein is based on polynomial chaos (PC) expansions, which have been widely used in stochastic mechanics since the development of the so-called stochastic spectral methods [24]. These approaches rely upon a Fourier–Hermite expansion of the uncertain input parameters and the model output into the governing equations. The coefficients of the PC representation of the model response are eventually determined using a Galerkin minimization technique. Such methods are commonly named intrusive since the solving scheme of the deterministic problem has to be modified. As an alternative, non-intrusive schemes have recently emerged in the literature, allowing the computation of the PC coefficients by means of a set of deterministic model evaluations, i.e. without modifying the underlying computer code. Two approaches are usually distinguished:

  • the projection approach: each PC coefficient is recast as a multidimensional integral [25], [26], [27], [28], [29], [30], [31] which can be computed either by simulation or quadrature methods;

  • the regression approach [29], [32], [33], [34]: the PC coefficients are estimated by minimizing the mean square error of the response approximation in the mean square sense.

The regression scheme reveals efficient when dealing with a moderate number of input variables [29], [32], [33], [34], [35]. It was successfully applied to global sensitivity analysis in [34] where the link between the orthogonal PC basis and the Sobol’ decomposition [7] was highlighted. However, one limitation of this technique lies in the polynomial increase of the required number of model evaluations (i.e. the computational cost) with the number of terms in the PC expansion, which itself dramatically increases with the number of input parameters. To circumvent this problem, two approaches based on sparse PC expansions have been recently proposed, namely a stepwise regression technique [36], [37] and a least angle regression algorithm [38]. Both strategies allow one to adaptively retain basis functions in the PC expansion.

The aim of the present paper is to show how the sparse PC expansions may be used in global sensitivity analysis and in particular illustrate their accuracy and great computational efficiency. In this respect, various analytical models will be addressed as well as a structural mechanics problem.

The remainder of this paper is organized as follows. In the next section one presents the framework of the ANOVA decomposition of the model response. Section 3 is dedicated to the PC representation of the model response and the estimation of the PC coefficients. Two estimates of the accuracy of the PC approximations are proposed and the truncation strategy of the representation is discussed. The iterative procedure for building up a sparse PC approximation is described in Section 4. PC-based sensitivity indices are proposed in Section 5. The computational gain provided by the sparse PC expansions compared to crude Monte Carlo simulation is illustrated in Section 6 by numerical examples, i.e. the studies of three analytical functions and the structural mechanics finite element model of a frame structure.

Section snippets

Probabilistic formulation of the problem

Let us consider a physical model represented by a deterministic function y=M(x). Here x={x1,,xM}TRM, M1 is the vector of the input variables, and y={y1,,yQ}TRQ, Q1 is the vector of quantities of interest provided by the model, referred to as the model response in the sequel. As the input vector x is assumed to be affected by uncertainty, it is represented by a random vector X with prescribed joint probability density function (PDF) fX(x). This leads to introduce the probability space (RM,B

Polynomial chaos representation

Let us first denote by L2i the Hilbert space L2(R,fXi) of square-integrable functions with respect to the probability measure fXi(xi)dxi, for all i=1,…,M. Let (πj(i))jN be a complete orthonormal basis (CONB) of L2i. The orthonormality property meansπj(i),πk(i)Li2E[πj(i)(Xi)πk(i)(Xi)]=δj,kwhere δj,k=1 if j=k and 0 otherwise. In this work a focus is given to polynomial bases. In practice, most common continuous distributions can be associated to a specific family of polynomials [44], [45],

Sparse polynomial chaos approximation

Let A be a non-empty finite subset of NM. It is possible to consider the truncated PC expansion defined byYMA(X)αAaαψα(X)The set A is referred to as the truncation set in the sequel. Note that the common truncation scheme in Eq. (16) corresponds to the choice A=AM,p, where p is a fixed integer. However, it has been shown that the high cardinality of this set (say card(A)=(M+pp)) may lead to an unaffordable computational cost.

Of interest is the determination of truncation sets A of low

ANOVA decomposition of the polynomial chaos expansion

Consider a sparse PC expansion of the model response produced by the procedure outlined in the last section:Y=M(X)MA(X)=αAa^αψα(X)The statistical moments of the response PC expansion can be analytically derived from its coefficients. In particular, the mean and the variance, respectively, readμA,Y=a^0DA=αA{0}a^α2Let us define by Ii1,,is the set of α-tuples in A such that only the indices {i1,…,is} are non-zero:Ii1,,is={αA:αk=0k(i1,,is),k=1,,M}The set Ii1,,is corresponds to the

Basic results

Let us consider the so-called Ishigami function which is widely used for benchmarking in global sensitivity analysis [1], [67]:Y=sinX1+7sin2X2+0.1X34sinX1where the Xi (i=1,…,3) are independent random variables that are uniformly distributed over [π,π]. Note that the model under consideration is sparse in nature since:

  • the maximal interaction order is 2;

  • the underlying function is even with respect to the variables X2 and X3, hence the PC terms associated with odd polynomials in these variables

Conclusion

In this paper, a method is proposed to compute sensitivity indices of the response of mathematical models at a low computational cost. It relies upon the substitution of the model under consideration by its polynomial chaos approximation. In order to reduce the number of unknown PC coefficients to identify and hence the required number of computer experiments (i.e. the computational cost), an adaptive algorithm is proposed for automatically detecting the significant PC terms, leading to a

References (75)

  • O.P. Le Maître et al.

    A stochastic projection method for fluid flow—II. Random process

    J Comput Phys

    (2002)
  • H.G. Matthies et al.

    Galerkin methods for linear and nonlinear elliptic stochastic partial differential equations

    Comput Methods Appl Mech Eng

    (2005)
  • P.S. Beran et al.

    Uncertainty quantification of limit-cycle oscillations

    J Comput Phys

    (2006)
  • B. Sudret

    Global sensitivity analysis using polynomial chaos expansions

    Reliab Eng Syst Safety

    (2008)
  • G. Blatman et al.

    Sparse polynomial chaos expansions and adaptive stochastic finite elements using a regression approach

    C R Mécanique

    (2008)
  • G. Blatman et al.

    An adaptive algorithm to build up sparse polynomial chaos expansions for stochastic finite element analysis

    Probab Eng Mech

    (2010)
  • T. Homma et al.

    Importance measures in global sensitivity analysis of nonlinear models

    Reliab Eng Syst Safety

    (1996)
  • A. Saltelli

    Making best use of model evaluations to compute sensitivity indices

    Comput Phys Commun

    (2002)
  • J.A.S. Witteveen et al.

    Modeling physical uncertainties in dynamic stall induced fluid–structure interaction of turbine blades using arbitrary polynomial chaos

    Comput Struct

    (2007)
  • J. An et al.

    Quasi-regression

    J Complexity

    (2001)
  • B. Ganapathysubramanian et al.

    Sparse grid collocation schemes for stochastic natural convection problems

    J Comput Phys

    (2007)
  • C. Tong

    Refinement strategies for stratified sampling methods

    Reliab Eng Syst Safety

    (2006)
  • C.J. Sallaberry et al.

    Extension of Latin hypercube samples with correlated variables

    Reliab Eng Syst Safety

    (2008)
  • I.M. Sobol’

    Theorems and examples on high dimensional model representation

    Reliab Eng Syst Safety

    (2003)
  • J.C. Helton et al.

    Computational implementation of sampling-based approaches to the calculation of expected dose in performance assessments for the proposed high-level radioactive waste repository at Yucca Mountain, Nevada

    Reliab Eng Syst Safety

    (2009)
  • P.-L. Liu et al.

    Optimization algorithms for structural reliability

    Struct Safety

    (1991)
  • D. Wei et al.

    Structural reliability analysis by univariate decomposition and numerical integration

    Probab Eng Mech

    (2007)
  • C. Xu et al.

    Uncertainty and sensitivity analysis for models with correlated parameters

    Reliab Eng Syst Safety

    (2008)
  • Saltelli A, Chan K, Scott EM. Sensitivity analysis. John Wiley & Sons;...
  • B. Efron et al.

    The Jacknife estimate of variance

    Ann Stat

    (1981)
  • H. Cukier et al.

    Nonlinear sensitivity analysis of multiparameter model systems

    J Comput Phys

    (1978)
  • A. Saltelli et al.

    A quantitative, model independent method for global sensitivity analysis of model output

    Technometrics

    (1999)
  • I.M. Sobol’

    Sensitivity estimates for nonlinear mathematical models

    Math Model Comput Exp

    (1993)
  • G.E. Archer et al.

    Sensitivity measures, ANOVA-like techniques and the use of bootstrap

    J. Stat. Comput. Simul.

    (1997)
  • A. Saltelli et al.

    A quantitative model-independent method for global sensitivity analysis of model output

    Technometrics

    (1999)
  • I.M. Sobol’ et al.

    Global sensitivity indices for nonlinear mathematical models. Review

    Wilmott Mag

    (2005)
  • M.D. Morris

    Factorial sampling plans for preliminary computational experiments

    Technometrics

    (1991)
  • Cited by (0)

    View full text