An adaptive high-dimensional stochastic model representation technique for the solution of stochastic partial differential equations

https://doi.org/10.1016/j.jcp.2010.01.033Get rights and content

Abstract

A computational methodology is developed to address the solution of high-dimensional stochastic problems. It utilizes high-dimensional model representation (HDMR) technique in the stochastic space to represent the model output as a finite hierarchical correlated function expansion in terms of the stochastic inputs starting from lower-order to higher-order component functions. HDMR is efficient at capturing the high-dimensional input–output relationship such that the behavior for many physical systems can be modeled to good accuracy only by the first few lower-order terms. An adaptive version of HDMR is also developed to automatically detect the important dimensions and construct higher-order terms using only the important dimensions. The newly developed adaptive sparse grid collocation (ASGC) method is incorporated into HDMR to solve the resulting sub-problems. By integrating HDMR and ASGC, it is computationally possible to construct a low-dimensional stochastic reduced-order model of the high-dimensional stochastic problem and easily perform various statistic analysis on the output. Several numerical examples involving elementary mathematical functions and fluid mechanics problems are considered to illustrate the proposed method. The cases examined show that the method provides accurate results for stochastic dimensionality as high as 500 even with large-input variability. The efficiency of the proposed method is examined by comparing with Monte Carlo (MC) simulation.

Introduction

Over the past few decades there has been considerable interest among the scientific community in studying physical processes with stochastic inputs. These stochastic input conditions arise from uncertainties in boundary and initial conditions as well as from inherent random material heterogeneities. Material heterogeneities are usually difficult to quantify since it is physically impossible to know the exact property at every point in the domain. In most cases, only a few statistical descriptors of the property variation or the property variation in small test regions can be experimentally determined. This limited information necessitates viewing the property variation as a random field that satisfies certain statistical properties/correlations. Several techniques, such as Karhunen–Loève (K-L) expansion [1] and model reduction methods [2], [3] have been proposed to generate these stochastic input models. This naturally results in describing the physical phenomena using stochastic partial differential equations (SPDEs).

In the past decade, there has been tremendous progress in posing and solving SPDEs with the methods used usually classified into three major groups. The first group refers to sampling methods. The most traditional one is the Monte Carlo (MC) method. Its convergence rate does not depend on the number of independent input random variables. Furthermore, MC methods are very easy to implement given a working deterministic code. However, the number of realizations required to acquire good statistics is usually quite large. The second group of methods consists of moment/perturbation methods, e.g. KL-based moment-equation approach [4], [5], [6]. These methods can deal with large number of inputs. However, they are limited to small fluctuations and do not provide high-order statistics of the solution. Significant emphasis is given recently on the non-perturbative methods. The first approach in this group for quantifying uncertainty is the spectral stochastic finite element method (SSFEM) [1] based on generalized polynomial chaos expansions (gPC) [7]. In this method, we project the dependent variables of the model onto a stochastic space spanned by a set of complete orthogonal polynomials and then a Galerkin projection scheme is used to transform the original stochastic problem into a system of coupled deterministic equations which can be solved by the finite element method. The gPC was successfully applied to model uncertainty propagation in random heterogenous media [8], [9], [10], [11], [12]. Error bounds and convergence studies [13] have shown that these methods exhibit fast convergence rates with increasing orders of expansions. Although this method can deal with large variations of the property, the coupled nature of the resulting equations for the unknown coefficients in the spectral expansion makes the solution of the stochastic problem extremely complex as the number of stochastic dimensions and/or the number of expansion terms increase, the so called curse of dimensionality. In fact, computational complexity of the problem increases combinatorially with the number of stochastic dimensions and the number of expansion terms. In addition, it is required to develop a stochastic simulator, which is a non-trivial task especially if the underlying PDEs have complicated nonlinear terms. Therefore, this method is often limited to a small number of inputs (1–10) .

There have been recent efforts to couple the fast convergence of the Galerkin methods with the decoupled nature of MC sampling, the so called stochastic collocation method based on sparse grids [14], [15], [16], [17], [18], [19]. Hereafter, this method is referred as conventional sparse grid collocation (CSGC) method. This framework represents the stochastic solution as a polynomial approximation. This interpolant is constructed via independent function calls to the deterministic problem at different interpolation points which are selected based on the Smolyak algorithm [20]. This method is also applied to model processes in random heterogenous media [2], [3], [21], [22]. However, there are several disadvantages. As is well known, the global polynomial interpolation cannot resolve local discontinuity in the stochastic space. Its convergence rate still exhibits a logarithmic dependence on the dimension. For high-dimensional problems, a higher-interpolation level is required to achieve a satisfactory accuracy. However, at the same time, the number of collocation points required increases exponentially for high-dimensional problems (>10) as shown in Fig. 1. Therefore, its computational cost becomes quickly intractable. This method is still limited to a moderate number of random variables (5–15). To this end, Ma and Zabaras [23] extended this methodology to adaptive sparse grid collocation (ASGC). This method utilizes local linear interpolation and uses the magnitude of the hierarchical surplus as an error indicator to detect the non-smooth region in the stochastic space and thus place automatically more points around this region. This approach results in further computational gains and guarantees that a user-defined error threshold is met. In [23], it is shown that the ASGC can successfully resolve stochastic discontinuity problems and solve stochastic elliptic problems up to 100 dimensions when the weights of each dimension are highly anisotropic and thus when the ASGC places more points only along the first few important dimensions. However, it is also shown in the paper that when the importance of each dimension weighs equally, ASGC cannot solve the problem accurately even with a moderate (21) stochastic dimensionality. In this case, the effect of ASGC is nearly the same as of the CSGC and thus the convergence rate deteriorates. As is well known, in realistic random heterogeneous media often we deal with a very small correlation length and this results in a rather high-dimensional stochastic space with nearly the same weights along each dimension. In this case, all the previously mentioned stochastic methods are obviously not applicable.

These modeling issues for high-dimensional stochastic problems motivate the use of the so called “High-Dimensional Model Representation” (HDMR) technique. Hereafter, this method is referred as conventional HDMR. It is a general set of quantitative assessment and analysis tools for capturing the high-dimensional relationships between sets of input and output model variables. It was originally developed as a methodology to create an efficient fully equivalent operational model of the original chemical systems [24], [25], [26], [27], [28], [29]. It expresses the model output as an additive hierarchical superposition of correlated functions with increasing numbers of input variables, i.e. 1, 2,  ,  up to the total number of input variables. A systematic mapping procedure between the inputs and outputs is prescribed to reveal the hierarchy of correlations among the input variables. At each new level of HDMR, higher-order correlated effects of the input variables are introduced. If these higher-order correlated effects between the input variables have negligible effect upon the output, the HDMR approximation is accurate enough by using only lower-order (usually up to the third-order) component functions. This is the case for most realistic physical systems and is the ansatz that the HDMR is based upon. Depending on the way that one determines the hierarchical component functions in HDMR, there are particularly two types of HDMR: ANOVA-HDMR and CUT-HDMR [24]. ANOVA-HDMR is the same as the analysis of variance (ANOVA) decomposition used in statistics and is useful for measuring the contributions of variance of each component function to the overall variance and therefore for computing the global sensitivity [30]. It involves high-dimensional integration and thus is computationally expensive. On the other hand, the CUT-HDMR expansion is a finite exact representation of the model output in the hyperplane passing through a reference point in the input variable space, which is different from gPC that has infinite terms. It only involves function evaluations at the sample points and is more computationally efficient than ANOVA-HDMR. Therefore, it is our focus in this work.

The model output can be considered as a function taking value over the parameter space supported by the input variables. Therefore, it is also a multivariate function approximation/interpolation method as well as a means to analyze the relevant statistics of the random output, which is the same idea as stochastic collocation method where the solution to the SPDEs is regarded as a stochastic function in the stochastic input space. It is thus straightforward to apply CUT-HDMR in the random space to construct the stochastic input–output mapping. The most important aspect of this method is the capability to find a way to numerically represent each component function of HDMR. In Ref. [31], CUT-HDMR is applied to a transport model to represent the stochastic response and MC analysis is used on the obtained approximation to obtain statistics. However, each CUT-HDMR component function was numerically represented as a low-dimensional look-up table over its variables. To obtain an approximate value, one needs to search and interpolate in the table. In [32], [33], CUT-HDMR is derived from a Taylor expansion and is used to find moments of the solution using the so called moment-based quadrature rule in stochastic mechanics. Although the name “dimension reduction method” is used in these papers, it is another form of HDMR. Later, the same author applied this method to fracture and reliability analysis, where component functions are interpolated using tensor product Lagrange polynomials and again MC analysis is used on the obtained approximation to obtain statistics [34], [35]. The author in [36] also applied CUT-HDMR to reliability analysis however the interpolant is constructed through moving least squares (MLS) approximation which is only limited up to three dimensions. Most of the above applications approximate each component function on a tensor-product uniform sampling space thus they are very expensive and not very accurate implementations. In addition, all the previous applications are still limited to relatively low-number of stochastic dimensions (<12).

In Ref. [37], the authors have applied CUT-HDMR, under the name “Anchored-ANOVA”, to find the mean and variance of the solution using multi-element probabilistic collocation method and integrating the resultant HDMR expansion term by term. In Ref. [38], the authors have applied the same method for the computation of high-dimensional integrals using quadrature schemes with applications to finance. They also developed a dimension-adaptive version of HDMR to find the important component functions and related it with the anisotropic sparse grid method. Motivated by their work, we develop a general framework to combine the strength from both HDMR and ASGC and apply it to uncertainty quantification. We first redefine the way to compute the error indicator in the formulation of ASGC. The use of the magnitude of the hierarchical surplus as the error indicator is too sharp and may result in non-terminating algorithms. The new error indicator incorporates the information from both the basis function and the surplus. This guarantees that the refinement will stop at a sufficient interpolation level. Then HDMR is used to decompose the high-dimensional stochastic problem into several lower-dimensional sub-problems. Each low-dimensional sub-problem is solved by ASGC in a locally-adaptive way within only related dimensions. In this way, an efficient low-dimensional stochastic reduced-order model is constructed and any model output in the stochastic space can be interpolated. By using ASGC, the interpolation of component functions is done efficiently by summing the corresponding hierarchical surplus and basis function as compared with that of using a numerical table. Mean and variance can also be obtained analytically through integrating the basis functions without any MC analysis. In addition, the ASGC provides a linear combination of tensor products chosen in such a way that the interpolation error is nearly the same as for full-tensor product (numerical table) in higher dimensions. In practice, HDMR is often truncated into a lower-order representation. However, for a very large stochastic dimension (>100), even second-order expansion will have too many component functions. For example, 125,251 component functions are needed for a second-order expansion of a 500-dimensional problem. Therefore, we need to find a way to construct only the important component functions. Motivated by the work in [38], in this paper, we also develop a dimension-adaptive version of HDMR to detect the important component functions. This is defined in two steps: First, the first-order HDMR is constructed and a weight associated with each term is defined to identify the most important dimensions. Then, higher-order components functions are constructed which consist only of these important dimensions. This method to our knowledge is the first approach which can solve high-dimensional stochastic problems by reducing the dimensions from truncation of HDMR and resolve low-regularity by local adaptivity through ASGC. It is noted here that the adaptivity in our paper is different from that used in [38]. The definition of the error indicator is not the same. The authors in [38] did not identify the important dimensions and thus there is a significant computational overhead in finding the important component functions for high-dimensional problems. It is also noted that in [33], [37], the CUT-HDMR is written in an explicit form which is not suitable for adaptive construction. On the other hand, in our paper, a recursive form of an HDMR component function is introduced which is the basis for our adaptive implementation.

This paper is organized as follows: In the next section, the mathematical framework of stochastic PDEs is formulated. In Section 3, the ASGC method for solving PDEs is briefly reviewed. In Section 4, we use HDMR to introduce a new method for solving SPDEs. The numerical examples are given in Section 5. Finally, concluding remarks are provided in Section 6.

Section snippets

Problem definition

In this section, we follow the notation in [23]. Let us define a complete probability space (Ω,F,P) with sample space Ω which corresponds to the outcomes of some experiments, F2Ω is the σ-algebra of subsets in Ω and P:F[0,1] is the probability measure. Also, let us define D as a d-dimensional bounded domain DRd(d=1,2,3) with boundary D. We are interested in finding a stochastic function u:Ω×DR such that for P-almost everywhere (a.e.) ωΩ, the following holds:L(x,ω;u)=g(x,ω),xD,andB(x;u)=h

Adaptive sparse grid collocation method (ASGC)

In this section, we briefly review the development of the ASGC strategy. For more details, the interested reader is referred to [23].

The basic idea of this method is to have a finite element approximation for the spatial domain and approximate the multi-dimensional stochastic space Γ using interpolating functions on a set of collocation points {Yi}i=1kΓ. Suppose we can find a finite element approximate solution u to the deterministic solution of the problem in Eq. (9), we are then interested

High-dimensional model representations (HDMR)

In this section, the basic concepts of HDMR are introduced following closely the notation in [24], [25], [29], [38]. For a detailed description of the theory applied to deterministic systems, the interesting reader may refer to [25].

Let f(Y):RNR be a real valued smooth multivariate stochastic function. Here, it is noted that f(Y) may be also a function of physical coordinate, f(Y,x). From now on, we will omit x to simplify the notation. HDMR represents f(Y) as a finite hierarchical correlated

Numerical examples

In this section, two sets of numerical examples are considered to illustrate the proposed HDMR technique. In the investigations below we consider both adaptive HDMR and conventional HDMR where the adaptivity here refers only to the truncation of the expansion. In all examples, the ASGC method is used for the calculation of the component functions. The ASGC method used alone refers to the adaptive sparse grid collocation in [23] applied directly (without HDMR expansion) to the interpolation of

Conclusions

In this paper, a novel adaptive high-dimensional stochastic model representation technique for solving high-dimensional SPDEs is introduced. This method applies HDMR in the stochastic space and decomposes the original N-dimensional problem into several low-dimensional sub-problems. This has been shown to be more efficient than solving the N-dimensional problem directly. Each sub-problem is solved using ASGC, where a proper error indicator is introduced to adaptively refine locally the

Acknowledgments

This research was supported by the Computational Mathematics Program of AFOSR (Grant F49620-00-1-0373), the Computational Mathematics Program of the NSF (Award DMS-0809062) and an OSD/AFOSR MURI–09 award on uncertainty quantification.

References (43)

Cited by (199)

View all citing articles on Scopus
View full text