Sparse polynomial chaos expansion based on Bregman-iterative greedy coordinate descent for global sensitivity analysis

Highlights

• A novel algorithm is proposed to build sparse PCE for global sensitivity analysis.

• The accuracy and efficiency are assessed on benchmarks and engineering examples.

• A detailed comparison is made with OMP, LAR and two adaptive algorithms.

• The algorithm provides a better tradeoff among accuracy, complexity and efficiency.

Abstract

Polynomial chaos expansion (PCE) is widely used in a variety of engineering fields for uncertainty and sensitivity analyses. The computational cost of full PCE is unaffordable due to the ‘curse of dimensionality’ of the expansion coefficients. In this paper, a novel methodology for developing sparse PCE is proposed by making use of the efficiency of greedy coordinate descent (GCD) in sparsity exploitation and the capability of Bregman iteration in accuracy enhancement. By minimizing an objective function composed of the ${\mathcal{l}}_{\text{1}}$ norm (sparsity) of the polynomial chaos (PC) coefficients and regularized ${\mathcal{l}}_{\text{2}}$ norm of the approximation fitness, the proposed algorithm screens the significant basis polynomials and builds an optimal sparse PCE with model evaluations much fewer than unknown coefficients. To validate the effectiveness of the developed algorithm, several benchmark examples are investigated for global sensitivity analysis (GSA). A detailed comparison is made with the well-established orthogonal matching pursuit (OMP), least angle regression (LAR) and two adaptive algorithms. Results show that the proposed method is superior to the benchmark methods in terms of accuracy while maintaining a better balance among accuracy, complexity and computational efficiency.

Introduction

Along with the ever-increasing complexity of computer models for engineering simulations, the inherent uncertainties of input data and model parameters are evolving rapidly. In this context, characterizing the uncertainties within a computer model is of great importance, and this has motivated the development of a variety of numerical techniques for the emerging field of uncertainty quantification. To efficiently quantify the effect of variation in each input parameter on model outputs, one popular technique is to substitute a computationally expensive model with a surrogate model that possesses similar quantities of interest such as statistical moments and the distribution of model outputs.

Surrogate model (also known as metamodel) is a mathematical or numerical approximation of a complex model generated by mapping from a small amount of random inputs to the corresponding model outputs. Over the past few years, a number of surrogate models have been developed in the field of uncertainty quantification, for example polynomial regression model [1], radial basis function [2], Kriging [3]/Gaussian process [4], artificial neural network [5], support vector regression [6], [7], ensemble of surrogates [8] and PCE [9], [10], [11], [12], [13], [14], [15], among which PCE has received much attention for uncertainty and sensitivity analyses [9], [10], [11], [12], [13], [14], [15], [16].

First introduced by Ghanem and Spanos [17] to stochastic mechanics based on homogeneous chaos theory [18] and later generalized by Xiu and Karniadakis [19] for different types of statistical distributions (e.g. uniform, beta and gamma), the PCE approach is to represent explicitly the stochastic model response as a series of orthonormal multivariate polynomials, i.e. PC basis [17]. In this scenario, quantification of the response probability density function is equivalent to estimation of the PC coefficients that are the coordinates of the stochastic response in the basis and can be evaluated at a set of sampling points in the input space. To build a PCE, two main approaches are typically adopted: intrusive and non-intrusive. The intrusive approach requires modifying the solving scheme of the deterministic governing equations of the model [20], whereas the non-intrusive approaches such as the projection method [21] and the regression method [9], [22] compute the PC coefficients by performing repeated simulations on limited number of input–output samples. Nevertheless, for either intrusive or non-intrusive approaches, the number of model evaluations (i.e. the computational cost) required for computing the PC coefficients increases dramatically with the number of input variables and the order of expanded polynomials.

To circumvent the issue of ‘curse of dimensionality’, several efficient non-intrusive approaches have been developed in recent years, such as adaptive methods [10], [23], [24], [25] for sequentially selecting the significant basis polynomials from the full PCE, multi-fidelity methods [4], [26] for achieving accurate predictions of quantities of interest using combination of ‘low-fidelity’ and ‘high-fidelity’ PCE simulations, sparse grid based methods [27], [28] for utilizing sparse grid interpolation techniques to reduce the number of collocation points in constructing the PCE, ${\mathcal{l}}_1$-minimization methods [26], [29], [30], [31] for scanning important bases by minimizing the ${\mathcal{l}}_1$ norm of PCE coefficients while preserving the fitting accuracy, and Bayesian methods [32], [33] for generating the sparse PCE with statistical model selection criteria. These approaches have been demonstrated with considerable computational gains compared to the classical PCE method.

In the field of signal processing and data analysis, coordinate descent (CD) algorithms have received much attention due to its simplicity and efficiency in sparsity exploitation [34], [35], [36], [37]. On the other hand, the least absolute shrinkage and selection operator (LASSO) is widely adopted to perform continuous model selection and enforce sparse solutions for problems where the number of predictors exceeds the number of cases [38]. It was found by Wu and Lange [35] that both cyclic CD, also known as pathwise coordinate descent (PCD) [34], and GCD were superior to the LAR [39] in terms of efficiency, robustness and model selection for LASSO-penalized ${\mathcal{l}}_2$ regression, while GCD was substantially faster than cyclic CD for LASSO-penalized ${\mathcal{l}}_1$ regression. GCD was further developed by Li and Osher [40] with combination of Bregman iteration [41] to solve the compressed sensing problems [42], [43], [44]. Recently, Zhou et al. [45] developed an adaptive method based on partial least squares and distance correlation for building sparse PCE and compared their algorithm with LASSO-based sparse PCE, in which PCD [36] was employed.

Motivated by the preceding analysis, the present paper aims at efficiently building sparse PCEs in the context of LASSO-based regression by taking advantage of both GCD and Bregman iteration. To the best of the authors’ knowledge, there is still no research on the capability of GCD algorithms for the purpose of sparse PCE construction. It has been shown that GCD may converge to a sparse solution significantly faster than cyclic or randomized CD [35], [46], especially for high-dimensional problems. In particular, GCD applied to problems of the LASSO form can sometimes approach to an optimal solution before executing even a single pass of all coordinates [46]. This suggests that GCD has an inherent screening ability for sparse optimization, which strongly motivates its combination with sparse PCE metamodeling for tackling high-dimensional problems. The novelty and contribution of the proposed method lie in the following aspects: (1) This study is probably the first work to develop a GCD algorithm for sparse PCE metamodeling in the context of LASSO-based regression. By updating the coordinate with the largest energy decrease [40], the GCD method is straightforward, efficient and robust in sparsity exploitation. (2) The regularization parameter $\lambda$ in the LASSO-based regression is commonly selected by cross-validation [35], which often leads to inefficient computation. In addition, the inefficiency is created by the dilemma that bigger $\lambda$ is preferred for more accuracy whereas smaller $\lambda$ gives rise to faster convergence in the coordinate updating of GCD [40]. To tackle this issue, this study incorporates Bregman iteration into the GCD to form a Bregman-iterative GCD (BGCD) algorithm structure, which not only settles the above inefficiency by adopting a moderate $\lambda$ with a relatively wide range of appropriate values, but also considerably enhances the convergence and accuracy of the GCD in solving the PC coefficients. (3) Existing non-intrusive algorithms for building sparse PCEs are mainly for problems with more cases than the number of predictors. However, the proposed algorithm is devised to handle problems with far fewer cases than the number of predictors, which is of great potential to the uncertainty quantification in practical engineering.

In this paper, the sparse PCE is employed for the GSA that aims at quantifying the respective effects of different inputs and their interactions on an assigned output response. GSA can provide complete information about the model behavior when the inputs vary in the entire domain. The priority level and ranking of the inputs resulting from the GSA can be very helpful for designers in narrowing the uncertain scope of model response. Among the developed works on sensitivity analysis, Sobol’ indices have attracted greater portion of attention due to the fact that they can provide accurate information for most models [16], [47], [48], [49]. The remainder of the paper is organized as follows. Section 2 gives a brief introduction of the Sobol’ decomposition and corresponding sensitivity indices for the GSA. In Section 3, the polynomial chaos approximation of a multidimensional model is recalled. The detailed procedure on how to determine the Sobol’ indices from the PC coefficients is given in particular. Then, in Section 4, a BGCD algorithm is proposed for building the optimal sparse PCE from a given sample set. The performance of the proposed algorithm is assessed on several benchmark examples for GSA in Section 5, and the results are compared with the well-established OMP algorithm, LAR technique and two adaptive methods. Section 6 summarizes concluding remarks.