Sparse polynomial surrogates for aerodynamic computations with random inputs
Abstract.
This paper deals with some of the methodologies used to construct polynomial surrogate models based on generalized polynomial chaos (gPC) expansions for applications to uncertainty quantification (UQ) in aerodynamic computations. A core ingredient in gPC expansions is the choice of a dedicated sampling strategy, so as to define the most significant scenarios to be considered for the construction of such metamodels. A desirable feature of the proposed rules shall be their ability to handle several random inputs simultaneously. Methods to identify the relative ”importance” of those variables or uncertain data shall be ideally considered as well. The present work is more particularly dedicated to the development of sampling strategies based on sparsity principles. Sparse multidimensional cubature rules based on general onedimensional GaussJacobitype quadratures are first addressed. These sets are non nested, but they are well adapted to the probability density functions with compact support for the random inputs considered in this study. On the other hand, observing that the aerodynamic quantities of interest (outputs) depend only weakly on the crossinteractions between the variable inputs, it is argued that only loworder polynomials shall significantly contribute to their surrogates. This ”sparsityofeffects” trend prompts the use of reconstruction techniques benefiting from the sparsity of the outputs, such as compressed sensing (CS). CS relies on the observation that one only needs a number of samples proportional to the compressed size of the outputs, rather than their uncompressed size, to construct reliable surrogate models. The results obtained with the test case considered in this work corroborate this expected feature.
Key words and phrases:
Computational fluid dynamics, Reynoldsaveraged NavierStokes equations, Uncertainty quantification, Polynomial chaos, Compressed sensingNomenclature
Polynomial basis  
Chord length  
Drag coefficient  
Lift coefficient  
Pitching moment coefficient  
Static pressure coefficient  
Total pressure coefficient  
Parameters space dimension  
Mathematical expectation  
Generic physical model  
Surrogate model  
Discretized computational model  
Polynomial chaos expansion coefficient of the computational model  
Thickness  
Freestream Mach number  
Measurement matrix  
Number of cubature points  
Level of the cubature rule  
Total order of the polynomial surrogates  
Number of polynomials  
Probability density function of the random input parameters  
Polynomial set of total order  
Thicknesstochord ratio  
Reynolds number  
Sparsity (the number of nonzero entries)  
Output quantities of interest  
Weighting matrix  
Angle of attack  
Beta law of the first kind  
Kurtosis  
Restricted isometry constant  
Polynomial truncation error  
Skewness  
Mean  
Mutual coherence  
th sensing vector  
Sensing basis  
th representation vector (polynomial chaos)  
Representation basis  
Standard deviation  
Onedimensional quadrature rule  
Cubature set  
Random input parameters  
th cubature weight  
th cubature point in the parameters space  
Subscript  
Index of the polynomial chaos  
Index of a cubature point 
1. Introduction
Surrogate models are typically used to perform optimization or uncertainty quantification (UQ) studies involving a complex modeling and simulation process, as encountered in computational fluid dynamics (CFD) among other engineering sciences applications. The principle of a surrogate model relies on an interpolation or regression procedure to estimate a scalar or vector field using a sampling data set constituted by carefully chosen outputs of the aforementioned complex process. Since the latter is often computationally expensive to run, the surrogate model shall be able to work out reliable estimates of its outputs at no extra computational costs except the evaluation of the surrogate itself. It should thus offer a non intrusive alternative to expensive runs of a complex process in order to sweep across the parameter space that influences it. When considering large parameter spaces, efficient algorithms are needed to provide an accurate surrogate representation of such parametric outputs. The kriging procedure [1] has gained a large attention over the past decades due to its robustness, accuracy, and ability to provide an estimate of the error done by the procedure; see the reviews by Kleijnen [2] or Bompard [3] and references therein. It has been applied to the simulation of air flow around a wing profile with the consideration of variable geometrical parameters in a previous study [4].
The polynomial chaos (PC) expansion is also a powerful tool for constructing spectrallike surrogate models of a parameterized process. A general methodology based on the Galerkin method has originally been proposed by Ghanem & Spanos [5] for the computation of the PC coefficients of the solution of a parameterized partial differential equation (PDE). This original scheme is heavily intrusive and has prompted the development of nonintrusive schemes, especially for PDEs arising in fluid dynamic models [6]. Two approaches for computing the coefficients of a PC expansion have usually been considered: (i) a projection approach, in which they are computed by structured (i.e. Gauss quadratures) or unstructured (i.e. MonteCarlo or quasi MonteCarlo sampling) quadratures; or (ii) a regression approach, minimizing some error measure or truncation tolerance. Both techniques suffer from some well identified shortcomings when the dimension of the parameter space, and the number of model evaluations alike, increase. Indeed, a PC expansion of total degree in variable parameters contains coefficients. A direct way to compute them is to use a tensor product grid in the parameter space requiring evaluations of the process, where the number of evaluations needed for one particular direction in that space is related to . These runs are very often unaffordable for large parameter spaces and complex configurations, as in CFD for example. Fortunately, the Smolyak algorithm [7] defines sparse grid quadratures involving points while preserving a satisfactory level of accuracy. Consequently, collocation techniques with sparse quadratures or adaptive regression strategies have been developed in order to circumvent this dimensionality concern [8, 9, 10, 11].
Other directions may be considered to deal with large parameter spaces. In the work presented in this paper we aim at benefiting from the sparsity of the process outputs themselves to reconstruct their PC representations in a nonadaptive way [12]. Indeed, we rely on the common observation that many crossinteractions between the input parameters are actually smoothened, or even negligible, once that have been propagated to some global quantities of interest processed from, say, complex aerodynamic computations. We can therefore expect to achieve a successful output recovery by the techniques known under the terminology of compressed sensing (CS) [13, 14]. In this theory the reconstruction of a sparse signal on a given, known basis requires only a limited number of evaluations at randomly selected points–at least significantly less than the a priori dimension of the basis. We thus resort to unstructured sampling sets to recover sparse outputs. We may also resort to highly structured sampling sets, such as nesting sets in large parameter spaces. The objective in this case is to be able to enrich the surrogate models by structured samples compatible with the structured samples used in a previous evaluation, without the need to reevaluate the whole process. Classical nested quadratures in one dimension go from a set of points to a set of points (for Fejér second rule [15]) or points (for ClenshawCurtis rule [16, 17]), such that nested rules involve sets of nodes of which dimensions double each time their socalled level is incremented, . These sets can subsequently by used with the Smolyak algorithm [7] to construct nested sparse sets in higher dimensions, together with adaptive strategies as developed by Gerstner & Griebel [18] for example.
We will not consider nested rules in this paper though, leaving these developments to another publication. As invoked above, we rather focus on sparse multidimensional cubature rules based on onedimensional Gausstype quadratures, on one hand. These sets are non nested, but they are well adapted to probability density functions with compact supports for the input parameters considered in the example addressed in this work. On the other hand, the ”sparsityofeffects” trend observed in complex simulations prompts the use of dedicated reconstruction techniques benefiting from the sparsity of the outputs. Since only loworder polynomials will significantly contribute to the output surrogates, one may infer from CS theory that only a number of samples proportional to the compressed size of the output, rather than the uncompressed size, is actually needed. In addition, these sampling points should be chosen randomly–i.e. they are unstructured. Therefore we may invoke sparsity patterns either for the input parameters or the output quantities of interest, to construct polynomial surrogates of the latter.
The rest of the paper is organized as follows. Section 2 introduces the basic configuration and CFD tools, namely a twodimensional RAE 2822 transonic airfoil [20, 21] and the Onera inhouse elsA software[22, 23]. This example serves as a guideline for the metamodeling techniques based on generalized PC expansions introduced in the subsequent sections: construction by quadrature rules (structured grids in the parameter space) in section 3 or by compressed sensing (using unstructured grids) in section 5, while the intermediate section 4 details the particular application of the former approach to the RAE 2822 airfoil. Some general conclusions and perspectives are drawn in section 6.
2. Problem setup
We start by introducing the problem setup and the numerical tools used to solve it. We consider a twodimensional transonic flow around an RAE 2822 transonic airfoil solved by steadystate ReynoldsAveraged NavierStokes (RANS) equations together with a SpalartAllmaras turbulence model closure [19]. The nominal flow conditions are the ones described in Cook et al. [20] for the test case #6 together with the wall interference correction formulas derived in [24, pp. 386–387] and their slight modifications suggested in [25, p. 130] (see also the CFD verification and validation website of the NPARC Alliance [21]). The nominal freestream Mach number , angle of attack , and Reynolds number (based on the airfoil chord length , fluid velocity, temperature and molecular viscosity at infinity) arise from the corrections and given in [25, p. 130] for the test case #6 outlined in Cook et al. [20], for which , , and .
2.1. Numerical nominal model
The nominal problem is discretized using a mesh shown in Fig. (1) and Fig. (2), where the boundary at infinity was left intensionally far (at about from the airfoil). These values proved to be sufficient to avoid spurious reflection with the farfield boundary. The discretized numerical solution is computed using the cellcentered finite volume CFD software elsA [22, 23]. It is based on:

First order Roe fluxes for the advection term of the turbulent variable;

Corrected second order diffusive terms based on the corrected mean of the cellcentered gradients of the two adjacent cells (referred to as the ”5p_cor” approach);

Source terms for the turbulent transport computed using the temperature gradients at the center of the cells.
The flow is attached with a weak shockwave on the suction side. The contour plot of the magnitude of the velocity are displayed on Fig. (3) and the static pressure profile at the wall are displayed on Fig. (4). Given the large number of simulations to run, the numerical parameters of the steady state algorithm proved to be essential to insure a fast convergence. This was performed using the following:

An implicite algorithm based on the LowerUpper Symmetric Successive Overrelaxation (LUSSOR) scheme [28] using relaxation cycles, increasing the CFL number after iterations to ;

A multigrid approach for the NavierStokes system over two grid levels with two iterations on the coarser grid;

A single fine level iteration for the turbulent equation alternating with a multigrid iteration for the RANS system.
Once the numerical parameters have been fixed, the number of iterations is determined from the evolution of the resulting global forces. A number of iterations (the discrete residuals of all equations and their decrease being checked at every iteration) appeared to be acceptable. Hence this number of iterations has been retained for all calculations so far. Further discussions on this issue are available in Dumont et al. [4] (available on demand).
2.2. Definition of the uncertainties
The aim of this work is to characterize the influence of uncertainties on the freestream Mach number , angle of attack , and thicknesstochord ratio on some aerodynamic quantities of interest, such as the drag, lift, or pitching moment coefficients , , or , respectively. These variable parameters are assumed to be independent and to follow Beta laws of the first kind . Therefore their probability density functions (PDF) read:
In the above is the usual Gamma function, and is the compact support of the random parameter . The parameters as well as the bounds for the three variable parameters , , are gathered in the Table 1 below. We note in passing that the model is the one arising from Jaynes’ maximum entropy principle [29, 30] when constraints on (i) the boundedness of the support, and (ii) the values of the averages and , are imposed. Here and in the subsequent developments the usual notation for is employed.
3. Polynomial surrogates
The computation of the first moments (mean, standard deviation, skewness, kurtosis…) of the aerodynamic quantities of interest when the variability of the parameters above is accounted for, is done thank to surrogate models, or response surfaces. We more particularly focus on polynomial surrogates in this study.
Let be a generic physical model involving parameters , such that the quantity of interest is given as:
Let be the set of dimensional polynomials with total order with respect to . We first note that this set has cardinality such that:
(1) 
A polynomial surrogate model of order for the model is obtained as:
(2) 
where is the marginal PDF of the random parameters with values in . The accuracy of this approximation may be assessed considering the limit of the meansquare norm as . However such a ”convergence” does not necessarily holds, and it depends on the probability measure .
There are otherwise several ways to construct polynomial surrogates: embedded projection (this is the original spectral stochastic finite element method of Ghanem & Spanos [5] which is highly intrusive), nonintrusive projection (or collocation) [10, 11], kriging [1, 2, 3], etc.; see also Le Maître & Knio [6] for a detailed introduction. Regression is also an alternative, whereby an optimization problem is formed. For that purpose, a set of sampling points is needed in order to discretize the minimization problem (2). It is assumed in the following that these points are chosen as the integration points of a cubature rule , which provides with positive weights and nodes in such that for a smooth function one can evaluate its average by:
(3) 
3.1. Regression approach
Using the foregoing cubature rule, the regression approach is formulated as a weighted leastsquares minimization problem for the coefficients of the polynomial surrogate of expanded on the monomials of partial total order (i.e. with ). Let where is given by Eq. (1), then:
where , , and ; also stands for the transpose of . Numerous methods are available to solve this problem whenever ; we do not follow this approach in the subsequent developments though. We are rather interested in the situation whereby , and more interestingly . It is addressed subsequently in the section 5.
3.2. Projection approach
We assume now that a polynomial basis of is available. Then we construct the polynomial surrogate of by standard projection on the finite dimensional subspace of spanned by the truncated family of orthonormal polynomials up to the total order denoted by , where is again given by Eq. (1). The orthonormalization of this basis reads:
(4) 
Then where , . Using the cubature rule of Eq. (3), these expansion coefficients are approximated by:
This corresponds to the approximation . Such expansions are referred to as polynomial chaos (PC) expansions in the dedicated literature, provided that the variable parameters follow a multidimensional Gaussian (normal) distribution [5, 6]. They are otherwise called generalized polynomial chaos (gPC) expansions for other distributions [31, 32].
3.3. Application to Uncertainty Quantification (UQ)
Once the polynomial surrogate model has been derived, the first moments and/or cumulants of the quantity of interest can be computed using the cubature rule and evaluations of the physical model at these nodes. Indeed, for a regular function on one can estimate a mean output functional by:
The mean is obtained for , the variance is obtained for , the skewness for , the kurtosis for , etc. More generally, the th moment is obtained for , and may be used to compute the characteristic function :
where by the causality principle (or transport of PDFs) for one has:
Sensitivity indices may be computed alike. Denoting by the set of indices corresponding to the polynomials of depending only on the –th variable parameter , the maineffect gPCbased Sobol’ indices are given by (see e.g. Sudret [33]):
owing to the normalization condition (4). More generally, if is the set of indices corresponding to the polynomials of depending only on the parameters , the –fold joint sensitivity indices are:
In the subsequent application to the twodimensional configuration described in section 2, we will consider the maineffect sensitivity indices and the –fold joint sensitivity indices .
4. Application to the twodimensional RAE 2822 transonic airfoil
From the previous analysis, we see that the main ingredients requested for the construction of polynomial surrogates are the cubature rule and the truncated polynomial basis , for integration nodes and a multidimensional polynomial total order . In addition we have here for the parameter space dimension. Owing to the choices made for the variable parameters considered for this case (see Table 1), we have and . Therefore the integration points should be chosen from a GaussJacobi quadrature rule, and the polynomial basis should be constituted by the multidimensional Jacobi polynomials which are orthogonal with respect to the weight function .
4.1. Polynomial basis
The polynomial basis adapted to the parameters PDF is constituted by the threedimensional Jacobi polynomials , , such that . They are computed by:
where is the family of onedimensional orthonormal Jacobi polynomials with respect to the weight function .
In the present study the polynomial surrogates constructed for the evaluation of the drag, lift and pitching moment coefficients , and , respectively, are truncated up to the total order . Therefore multidimensional Jacobi polynomials are ultimately retained in those gPC expansions.
4.2. Cubature rules
Onedimensional GaussJacobi quadratures based on integration points are tailored to integrate on a smooth function :
(5) 
such that this rule turns to be exact for polynomials up to the order . Here is the number of fixed nodes of the rule, typically the bounds . Depending on the choice of , different terminologies are used:

is the classical GaussJacobi rule;

is the GaussJacobiRadau (GJR) rule, choosing or for instance;

is the GaussJacobiLobatto (GJL) rule, choosing and for instance.
Since in our case we have chosen a total order , GJL points are needed to recover exactly the orthonormality property (4) for the corresponding onedimensional Jacobi polynomials. Indeed should be defined such that in this situation. The –points GaussJacobi rules are illustrated on Fig. (5) for various values of the parameters of the Jacobi weight function , and the –points GaussJacobiLobatto rules are illustrated on Fig. (6). The blue dots correspond to and thus pertain to the PDF. The reason why we include the boundary nodes in the quadrature rule stems from the fact that the basic engineering practice would consider the evaluation of the physical model at the bounds of the support of the variable parameters. The main advantage of using GaussJacobi quadratures is that they do not add integration points for the increased order induced by the weight function . The ClenshawCurtis rule [16] for example is typically suited for polynomials of order , though higher orders are actually achieved in practice [17]. Thus if one uses nodes from this rule to compute the lefthand side of (5) an exact evaluation is achieved for polynomials up to the order , instead of with a GJL rule. However, the latter does not have the nesting property of the former.
Multidimensional quadratures (cubatures) may subsequently be obtained by tensorization and/or sparsification of these onedimensional rules. Firstly, a multidimensional tensorized grid is obtained by the straightforward product rule:
(6) 
It contains grid points, i.e. for the present Test Case (, ). Of course the order could be adapted depending on the –th variable parameter. Secondly, a sparse cubature rule can be derived thank to the Smolyak algorithm [7]. The socalled –th level, dimensional sparse rule is obtained by the following linear combination of product formulas [34]:
(7) 
For example, if and one has for the total number of grid points using a onedimensional GJL quadrature as the generating rule, and:
By a direct extension of the arguments devised by Novak & Ritter [35] or Heiss & Winschel [36], it can be shown that such a –th level, dimensional sparse rule based on the GJL onedimensional rule is exact for dimensional polynomials of total order . The total number of integration points of the rule is given for by the estimate (see e.g. Heiss & Winschel [36] and references therein), which typically outperforms the product rule with for . It is gathered in the Table 2 below for various combinations . The sparse rule is plotted in Fig. (8) for and , and in Fig. (10) for and together with the corresponding product rule in Fig. (7) and Fig. (9), respectively, for illustration purposes. We note that since the underlying onedimensional GJL rule is not nested, the multidimensional rules are not either. Also the weights of the latter may be negative for some nodes although the underlying onedimensional rules have positive weights.
\backslashbox  2  3  4  5  6 

2  4  8  16  32  64 
3  8  20  48  112  256 
4  17  50  136  352  880 
5  29  99  304  872  2384 
6  53  201  673  2082  6092 
7  85  363  1337  4483  14072 
8  133  647  2585  9293  31025 
9  193  1079  4697  18143  64469 
10  273  1769  8321  34323  129197 
4.3. Results
The first four moments of the drag, lift, and pitching moment coefficients , , and , respectively, are gathered in the Table 3 below for calls to the model using the elsA software [22, 23] with the –th level full product rule (6). Table 4 gathers the same results using a th level sparse rule (7) for which from Table 2. The reasons why we have not used a –th level sparse rule for a consistent comparison with the –th level product rule are twofold. Firstly, for this case , which is not competitive with the product rule. Secondly, the gPC coefficients are observed to be sparse so that the higherorder threedimensional polynomials contribute only marginally to the surrogates. This observation is elaborated further on in the subsequent section 5 dealing with the sparse reconstruction approach we have applied invoking the theory of compressed sensing. Using the th level sparse rule we are able to integrate exactly threedimensional polynomials up to total order , hence reconstruct polynomial surrogates up to total order , corresponding to multidimensional Jacobi polynomials in those gPC expansions.
133.37e04  34.128e04  1.0237e+00  3.3030e+00  
72.274e02  1.6695e02  9.6221e01  3.0630e+00  
453.99e04  32.239e04  5.7845e01  2.3190e+00 
133.38e04  34.097e04  1.0293e+00  3.2611e+00  
72.269e02  1.6729e02  9.8175e01  2.8678e+00  
453.96e04  32.175e04  5.9533e01  2.3885e+00 
The maineffect sensitivity indices computed with the –th level product rule are gathered in Table 5 below, while the joint sensitivity indices are gathered in Table 6. Tables 7 and 8 display the same indices computed with the –th level sparse rule. It is clearly apparent from these results that the freestream Mach number is the chief parameter controlling the variability of the aerodynamic coefficients in the range of analysis considered for this test case. We may also emphasize the discrepancies between the sensitivity indices computed with a full and a sparse rule, although their theoretical accuracy are different in the present case. These differences are even more pronounced for the joint sensitivity indices. The sparse reconstruction approach outlined in the next section yields sensitivities very comparable to the ones derived with the product rule, at a much lower computational cost though.
0.081e01  9.892e01  0.008e01  
0.034e01  9.554e01  0.286e01  
0.269e01  9.721e01  0.000e01 
0.021e02  0.000e02  0.172e02  
0.036e02  0.000e02  1.221e02  
0.007e02  0.000e02  0.089e02 
0.052e01  6.635e01  0.007e01  
0.033e01  5.883e01  0.195e01  
0.256e01  9.227e01  0.002e01 
0.138e02  23.944e02  0.847e02  
0.287e02  29.990e02  1.510e02  
0.085e02  3.276e02  0.582e02 
The PDFs of the three aerodynamic coefficients considered in this study are displayed on Fig. (11) through Fig. (13) and Fig. (14) through Fig. (16) using the –th level product rule and the –th level sparse rule, respectively. They were estimated from evaluations of the gPC surrogates and smoothing out the resulting histograms by a normal kernel density function. The means are shown on the plots with vertical blue lines.
Finally, the mean total pressure coefficients along the profile are displayed on Fig. (17) using the –th level product rule and on Fig. (18) using the –th level sparse rule. Error bars at , where is the standard deviation of the ’s, are also shown, together with the nominal total pressure coefficient (dotted line) obtained from a computation with the nominal parameters , , , and the nominal thicknesstochord ratio. We observe an unexpected drop of the standard deviation at locally at the suction side when the –th level sparse rule is used. This may be related to the negativeness of some weights with sparse rules, however we have not been able to find a more detailed explanation of this probable anomaly so far.
5. Nonadapted sparse reconstruction by –minimization
The results of the previous section, especially the maineffect and joint sensitivity indices, suggest that only loworder interactions exist between the variable parameters for their effect on the aerodynamic coefficients of interest. This indicates that among the gPC coefficients to be computed for the construction of their polynomial surrogates , the ones pertaining to the onedimensional polynomials (depending on a single variable parameter) dominate the others. Hence the vector of gPC coefficients of has many negligible components, so that it is compressible in the terminology of the theory of compressed sensing [13, 14, 37]. In this setting it is argued that only a number of samples proportional to the compressed size, rather than the uncompressed size, of the sought signal is needed in order to reconstruct it. This observation challenges the conventional approaches to sampling or imaging according to Shannon’s theorem, which states that the sampling rate (the Nyquist rate) must be at least twice the maximum frequency present in the signal. Compressed sensing, or compressive sampling (CS), asserts that one can recover some signals from far fewer samples or measurements than traditionally used in the widespread signal acquisition techniques. CS relies on two principles to make this possible:

Sparsity, which express the fact that many signals may have a concise representation once they are expressed in a proper basis ;

Incoherence, which express the fact that signals having a sparse representation in a given basis are actually spread out in the domain in which they are acquired. Or in other words, the sensing functions used to probe the signal have a dense representation in the basis .
The reconstruction procedure consists in correlating the signal with a small number of predefined sensing functions (for example, sinusoids if one aims at computing a Fourier transform) which are incoherent with the basis in which the signal is sparse. It is non adapted because it identifies the sparsity pattern, that is the order (location) of the negligible components of the signal in its sparsifying basis, and the leading components at the same time. The procedure can therefore efficiently capture the relevant information of a sparse signal without trying to comprehend that signal [37]. This is clearly a much desirable feature for practical industrial applications.
5.1. Theoretical background
We basically follow the introductory paper of Candès & Wakin [37] in this introductory section to CS. The sparse signal to be reconstructed is the polynomial surrogate model of total order , for which information is obtained by recording values of the model :
(8) 
Typically is a sinusoid and is a Fourier coefficient if a Fourier transform is applied, or a Dirac function if the model is evaluated at some sampling point : . This latter sensing procedure is the one considered in this work. Letting be the sensing matrix of which rows are the vectors ( is the conjugate transpose of ), the process of recovering a discretized version of the model from the observation of outputs reads:
This problem is generally illposed whenever , but CS theory tells us that a unique solution may be obtained if the vector is sparse, though. This is actually the case once the model is expanded on the orthonormal basis constituted by the polynomial chaos pertaining to the randomly variable parameters introduced in section 2. This property is observed a posteriori for the present problem from the results expounded in the foregoing section 4. Incidentally, it should be noted that sparsity can be proved a priori for some parametric, possibly nonlinear elliptic and parabolic partial differential equations in a general framework; see Chkifa et al. [38] and references therein. However, the extension of these analyses to the present RANS system supplemented with a turbulence closure model, not to say the NavierStokes equations, does not seem actually feasible.
The discretized version of the model in terms of the polynomial surrogate reads in its sparse representation on the basis . Applying the sensing procedure (8) to the polynomial surrogate yields:
(9) 
where is the sparse vector of the gPC coefficients in the basis to be computed, and is the socalled measurement matrix given by . Using as done here, it is a Vandermondetype matrix with . Again, Eq. (9) is illposed whenever , unless the sparsity of the sought solution is accounted for. Regularized versions of Eq. (9) exists for this case, which in turn ensure its wellposedness. Since the polynomial chaos expansion truncated at the total order is not complete for the exact representation of , a truncation error has also to be accounted for in the solution process. Together with the sparsity of the gPC coefficients, this can be accommodated by reformulating Eq. (9) as the following –minimization problem, known as Basis Pursuit Denoising (BPDN) [39]:
for some tolerance on the polynomial chaos truncation. The norms above are defined by , and is some diagonal weighting matrix. Its role is to prevent the algorithm from biasing toward the nonnegligible entries of of which associated columns in have large norms [12, 40]. Now the strategy for our present study is to solve with runs of the physical model significantly lower than the number of coefficients to be identified. CS shows that it is achievable provided that the target is actually sparse, or nearly sparse (compressive), and some constraints on the measurement matrix are fulfilled.
As already stated above, a key requirement for the successful recovery of a sparse signal is incoherence between the sensing basis and the representation basis . It is quantified by the following mutual coherence :
(10) 
where stands for the –th column of . It is a measure of how close to orthogonality the measurement matrix is. Based on this coherency measure, the following theorem from Candès & Romberg [41] asserts that if is sufficiently sparse in , the recovery of its gPC coefficients by –minimization is exact.
Theorem 5.1.
Assume that is –sparse on the gPC basis , that is it has at most nonzero entries. Then if sampling points are selected at random to form the measurement matrix , and:
for some constant , the solution of is exact with ”overwhelming” (sic) probability.
More precise results with structured random measurement matrices are given by e.g. Rauhut & Ward [45]. It should be noted that the role of coherence in this result is transparent. The smaller the coherence is, the closer the measurement matrix is to a unitary matrix, and the fewer sampling points are needed.
The previous Theorem 5.1 is however not entirely satisfactory from a practical point of view because (i) it does not allow for some truncation error, or noisy/imprecise measurements; (ii) it does not deal with approximately sparse (compressive) signals, for which a large subset of entries are negligible rather than strictly zero. These shortcomings may be alleviated simultaneously as established by Candès et al. [13]. To achieve this, a constraint on the measurement matrix needs be added to gain robustness in CS, the socalled restricted isometry property (RIP, also quoted as uniform uncertainty principle). For each integer , the isometry constant of is defined as the smallest number such that:
for all sparse vectors . Then is said to satisfy the RIP of order if, say, is not too close to . This property amounts to saying that all –column submatrices of are numerically wellconditioned, or columns selected arbitrarily in are nearly orthogonal. The following theorem by Candès et al. [13, 37] then states that can be solved efficiently:
Theorem 5.2.
This result calls for several comments. First, it is more general than Theorem 5.1 since, if the signal is exactly –sparse, and the reconstruction is exact whenever (noiseless case). Second, it deals with all signals, not the –sparse ones solely. Third, it is deterministic and does not involve any probability. Lastly, the bound on is the one originally proposed by Candès & Wakin [37] but it can be improved as proposed by e.g. Mo & Li [42]; such improvements are an active field of research at present.
5.2. Application to the twodimensional RAE 2822 transonic airfoil
We now apply the foregoing CS procedure to the nonadaptive computation of the gPC coefficients of the surrogate models for the aerodynamic coefficients , and . We use sampling points drawn at random following PDFs as defined in section 2. This sampling set is displayed on Fig. (19) below. The sole reason why we have chosen this sampling size if for its ease of use with the multithreading setup of our CFD software elsA [22, 23].
We subsequently apply BPDN to compute . For that purpose we use the Spectral Projected Gradient Algorithm (SPGL1) developed by van den Berg & Friedlander [43] and implemented in the Matlab package SPGL1 [44] to solve this –minimization problem. The tolerance was fixed at and we were able to find a solution for all surrogates with this a priori choice without resorting to crossvalidation, for example. Further investigations should be carried on on this topic, though. It should also be noted that no particular sampling strategy, such as stratification, lowdiscrepancy series, or preconditioning, has been applied at this stage to select the sampling set. Moreover, the weighting matrix was selected as the identity matrix. Alternative sampling and weighting strategies are outlined in several recent works[12, 40, 45, 46, 47, 48], though.
The mean and standard deviation of the drag, lift, and pitching moment coefficients , , and , respectively, are gathered in the Table 9 below. The maineffect sensitivity indices are gathered in Table 10, while the joint sensitivity indices are gathered in Table 11. These results are very close to the ones obtained by the –th level product rule. The PDFs of the aerodynamic coefficients considered here are displayed on Fig. (20) through Fig. (22). As for Fig. (11) through Fig. (16), they were estimated from evaluations of the gPC surrogates and smoothing out the resulting histograms by a normal kernel density function. The means are again shown on the plots with vertical blue lines. We finally compare on Fig. (23) through Fig. (25) (in scale) the PDFs computed by the three approaches considered in this work. We observe that the –th level product rule (with structured sampling points) and –minimization (with randomly selected sampling points) yield comparable results, but of course at a much lower computational cost with this latter technique.
133.33e04  34.052e04  
72.271e02  1.6703e02  
453.95e04  32.180e04 
0.080e01  9.893e01  0.008e01  
0.032e01  9.549e01  0.289e01  
0.268e01  9.722e01  0.000e01 
0.022e02  0.000e02  0.166e02  
0.031e02  0.003e02  1.265e02  
0.007e02  0.000e02  0.093e02 
6. Conclusions
In this paper we have addressed various methodologies with relevance to the construction of polynomial surrogates for parameterized complex processes as encountered in CFD problems. The present work has been more particularly focused on the development of dedicated sampling sets in the parameter space using either structured or unstructured grids. These techniques were illustrated with the example of an RAE 2822 airfoil in the transonic regime considering variable geometrical (the thicknesstochord ratio) and operational (the freestream Mach number and angle of attack) parameters.
Firstly, multidimensional sparse cubature rules based on onedimensional GaussJacobi rules have been used for uncertainty quantification of this twodimensional aerodynamic computation. The quantities of interest are the usual drag, lift, and pitching moment coefficients for which polynomial surrogates are sought for using the aforementioned sampling sets as learning sets. More particularly, GaussJacobiLobatto points have been considered since the probability density functions of the variable parameters have finite supports. Indeed, the engineering practice would typically include the boundary points of the parameter space in the learning sets.
Secondly, observing a posteriori that the aerodynamic quantities of interest are sparse in that parametric space, when projected on the multidimensional orthogonal polynomials associated to the parameters probability density functions, an –minimization procedure has been applied in the framework of the theory of compressed sensing. The latter allows to recover the expansion coefficients of the quantities of interest at a much lower computational cost than the sparse grids addressed in the first approach. Unstructured sampling points are needed in this process, selecting them randomly in the parameter space. Their number is typically less than the dimension of the polynomial space where the surrogates are sought for, and thus typically much less than the number of points of the multidimensional sparse rules that may be used for a given level of accuracy. The –minimization procedure is nonadaptive in the sense that it identifies both the amplitude of the leading expansion coefficients and their order. It thus constitutes a promising direction for future developments in practical applications for more complex geometries and flows, where adaptive strategies within the parametric space, weighted minimization, or preconditioned sampling sets may be needed.
Acknowledgments
This work has been supported by the European Union’s Seventh Framework Programme for research, technological development and demonstration under grant agreement #ACP3GA2013605036 (UMRIDA Project www.umrida.eu).
References
 Krige, D.G. A statistical approach to some basic mine valuation problems on the Witwatersrand. J. Chem. Metal. Mining Soc. South Africa, 52(6):119139 (1951).
 Kleijnen, J.P.J. Kriging metamodeling in simulation: A review. Eur. J. Operat. Research, 192(3):707716 (2009).
 Bompard, M. Modèles de substitution pour l’optimisation globale de forme en aérodynamique et méthode locale sans paramétrisation. PhD Thesis, Université NiceSophia Antipolis, December 2011.
 Dumont, A., Passaggia, P.Y., Peter, J., Salah el Din, I., Savin, É. Metamodels for the numerical errors with non intrusive UQ methods. UMRIDA Technical Report D2.312(9), Onera, Châtillon (2014).
 Ghanem, R., Spanos, P.D. Stochastic Finite Elements: A Spectral Approach. Springer, New York NY (1991).
 Le Maître, O., Knio, O. Spectral Methods for Uncertainty Quantification. With Applications to Computational Fluid Dynamics. Springer, Dordrecht (2010).
 Smolyak, S.A. Quadrature and interpolation formulas for tensor products of certain classes of functions. Soviet Math. Dokl., 4:240243 (1963).
 Eldred, M., Burkardt, J. Comparison of nonintrusive polynomial chaos and stochastic collocation methods for uncertainty quantification. AIAA paper #2009976 (2009).
 Blatman, G., Sudret, B. An adaptive algorithm to build up sparse polynomial chaos expansions for stochastic finite element analysis. Prob. Engng. Mech., 25(2):183197 (2010).
 Ghiocel, D.M., Ghanem, R.G. Stochastic finiteelement analysis of seismic soilstructure interaction. ASCE J. Engng. Mech., 128(1):6677 (2002).
 Xiu, D., Hesthaven, J.S. Highorder collocation methods for differential equations with random inputs. SIAM J. Sci. Comput., 27(3):11181139 (2005).
 Doostan, A., Owhadi, H. A nonadapted sparse approximation of PDEs with stochastic inputs. J. Comput. Phys., 230(8):30153034 (2011).
 Candès, E.J., Romberg, J.K., Tao, T. Stable signal recovery from incomplete and inaccurate measurements. Commun. Pure Appl. Math., 59(8):12071223 (2006).
 Donoho, D.L. Compressed sensing. IEEE Trans. Inform. Theory, 52(4):12891306 (2006).
 Fejér, L. Mechanische Quadraturen mit positiven Cotesschen Zahlen. Math. Z., 37(1):287309 (1933).
 Clenshaw, C.W., Curtis., A.R. A method for numerical integration on an automatic computer. Numer. Math., 2(1):197205 (1960).
 Trefethen, L.N. Is Gauss quadrature better than ClenshawCurtis? SIAM Rev., 50(1):6787 (2008).
 Gerstner, T., Griebel, M. Numerical integration using sparse grids. Numer. Algorithms, 18(34):209232 (1998).
 Spalart, P., Allmaras, S. A oneequation turbulence model for aerodynamic flows. AIAA Paper #19920439 (1992).
 Cook, P.H., McDonald, M.A., and Firmin, M.C.P. Aerofoil RAE 2822 – Pressure distributions, and boundary layer and wake measurements. In Experimental Data Base for Computer Program Assessment. AGARD Advisory Report No. 138, NATO, May 1979; Appendix A6.
 http://www.grc.nasa.gov/WWW/wind/valid/raetaf/raetaf01/raetaf01.html. Online the April 21, 2015.
 Cambier, L., Heib, S., Plot, S. The Onera elsA CFD software: input from research and feedback from industry. Mechanics & Industry, 14(3):159174 (2013).
 http://elsa.onera.fr/. Online the May 30, 2015.
 Garner, H.C., Rogers, E.W.E., Acum, W.E.A., and Maskell, E.C. Subsonic Wind Tunnel Wall Corrections. AGARDograph No. 109, NATO, October 1966; Chapter 6.
 Haase, W., Bradsma, F., Elsholz, E., Leschziner, M., and Schwamborn, D. (eds.). EUROVAL–An European initiative on validation of CFD codes. Notes on Numerical Fluid Mechanics 42, Vieweg Verlag, Wiesbaden (1993); Section 5.1.
 van Leer, B. Towards the ultimate conservative difference scheme, V. A second order sequel to Godunov’s method. J. Comput. Phys., 32(1):101136 (1979).
 van Albada, G.D., van Leer, B., Roberts, W.W. A comparative study of computational methods in cosmic gas dynamics. Astronomy & Astrophysics 108(1):7684 (1982).
 Yoon, S.K., Jameson, A. An LUSSOR scheme for the Euler and NavierStokes equations. AIAA Paper #19870600 (1987).
 Jaynes, E.T. Information theory and statistical mechanics. Phys. Rev., 106(4):620630 (1957).
 Jaynes, E.T. Information theory and statistical mechanics II. Phys. Rev., 108(2):171190 (1957).
 Xiu, D., Karniadakis, G.E. The Wiener–Askey polynomial chaos for stochastic differential equations. SIAM J. Sci. Comput., 24(2):619644 (2002).
 Soize, C., Ghanem, R. Physical systems with random uncertainties: chaos representations with arbitrary probability measure. SIAM J. Sci. Comput., 26(2):395410 (2004).
 Sudret, B. Global sensitivity analysis using polynomial chaos expansions. Reliab. Eng. Syst. Safe., 93(7):964979 (2008).
 Wasilkowski, G.W., Woźniakowski, H. Explicit cost bounds of algorithms for multivariate tensor product problems, J. Complexity, 11(1):156 (1995).
 Novak, E., Ritter, K. Simple cubature formulas with high polynomial exactness. Constructive Approx. 15(4):499522 (1999).
 Heiss, F., Winschel, V. Likelihood approximation by numerical integration on sparse grids. J. Econometrics, 144(1):6280 (2008).
 Candès, E.J., Wakin, M.B. An introduction to compressive sampling. IEEE Sig. Proc. Mag., 25(2):2130 (2008).
 Chkifa, A., Cohen, A., Schwab, C. Breaking the curse of dimensionality in sparse polynomial approximation of parametric PDEs. J. Math. Pure Appl., 103(2):400428 (2014).
 Chen, S.C., Donoho, D.L., Saunders, M.A. Atomic decomposition by basis pursuit. SIAM J. Sci. Comput., 20(1):3361 (1998).
 Peng, J., Hampton, J., Doostan, A. A weighted –minimization approach for sparse polynomial chaos expansions. J. Comput. Phys., 267:92111 (2014).
 Candès, E., Romberg, J.K. Sparsity and incoherence in compressive sampling. Inverse Probl., 23(3):969985 (2007).
 Mo, Q., Li, S. New bounds on the restricted isometry constant . Appl. Comput. Harm. Anal., 31(3):460468 (2011).
 van den Berg, E., Friedlander, M.P. Probing the Pareto frontier for basis pursuit solutions. SIAM J. Sci. Comput., 31(2):890912 (2008).
 van den Berg, E., Friedlander, M.P. SPGL1: A solver for largescale sparse reconstruction. June 2007.
 Rauhut, H., Ward, R. Sparse Legendre expansions via –minimization. J. Approx. Th., 164(5):517533 (2012).
 Yan, L., Guo, L., Xiu, D. Stochastic collocation algorithms using –minimization. Int. J. Uncertainty Quantification, 2(3):279293 (2012).
 Yuang, X., Karniadakis, G.E. Reweighted minimization method for stochastic elliptic differential equations. J. Comput. Phys., 248:87108 (2013).
 Hampton, J., Doostan, A. Compressive sampling of polynomial chaos expansions: Convergence analysis and sampling strategies. J. Comput. Phys., 280:363386 (2015).