Adaptive Smolyak Pseudospectral Approximations

Adaptive Smolyak Pseudospectral Approximations

Patrick R. Conrad and Youssef M. Marzouk

Polynomial approximations of computationally intensive models are central to uncertainty quantification. This paper describes an adaptive method for non-intrusive pseudospectral approximation, based on Smolyak’s algorithm with generalized sparse grids. We rigorously analyze and extend the non-adaptive method proposed in [6], and compare it to a common alternative approach for using sparse grids to construct polynomial approximations, direct quadrature. Analysis of direct quadrature shows that errors are an intrinsic property of some configurations of the method, as a consequence of internal aliasing. We provide precise conditions, based on the chosen polynomial basis and quadrature rules, under which this aliasing error occurs. We then establish theoretical results on the accuracy of Smolyak pseudospectral approximation, and show that the Smolyak approximation avoids internal aliasing and makes far more effective use of sparse function evaluations. These results are applicable to broad choices of quadrature rule and generalized sparse grids. Exploiting this flexibility, we introduce a greedy heuristic for adaptive refinement of the pseudospectral approximation. We numerically demonstrate convergence of the algorithm on the Genz test functions, and illustrate the accuracy and efficiency of the adaptive approach on a realistic chemical kinetics problem.

Key words. Smolyak algorithms, sparse grids, orthogonal polynomials, pseudospectral approximation, approximation theory, uncertainty quantification

AMS subject classifications. 41A10, 41A63, 47A80, 65D15, 65D32

1 Introduction

A central issue in the field of uncertainty quantification is understanding the response of a model to random inputs. When model evaluations are computationally intensive, techniques for approximating the model response in an efficient manner are essential. Approximations may be used to evaluate moments or the probability distribution of a model’s outputs, or to evaluate sensitivities of model outputs with respect to the inputs [21, 41, 33]. Approximations may also be viewed as surrogate models to be used in optimization [22] or inference [23], replacing the full model entirely.

Often one is faced with black box models that can only be evaluated at designated input points. We will focus on constructing multivariate polynomial approximations of the input-output relationship generated by such a model; these approximations offer fast convergence for smooth functions and are widely used. One common strategy for constructing a polynomial approximation is interpolation, where interpolants are conveniently represented in Lagrange form [1, 42]. Another strategy is projection, particularly orthogonal projection with respect to some inner product. The results of such a projection are conveniently represented with the corresponding family of orthogonal polynomials [4, 21, 43]. When the inner product is chosen according to the input probability measure, this construction is known as the (finite dimensional) polynomial chaos expansion (PCE) [14, 32, 8]. Interpolation and projection are closely linked, particularly when projection is computed via discrete model evaluations. Moreover, one can always realize a change of basis [9] for the polynomial resulting from either operation. Here we will favor orthogonal polynomial representations, as they are easy to manipulate and their coefficients have a useful interpretation in probabilistic settings.

This paper discusses adaptive Smolyak pseudospectral approximation, an accurate and computationally efficient approach to constructing multivariate polynomial chaos expansions. Pseudospectral methods allow the construction of polynomial approximations from point evaluations of a function [4, 3]. We combine these methods with Smolyak’s algorithm, a general strategy for sparse approximation of linear operators on tensor product spaces, which saves computational effort by weakening the assumed coupling between the input dimensions. Gerstner & Griebel [13] and Hegland [16] developed adaptive variants of Smolyak’s algorithm for numerical integration and illustrated the effectiveness of on-the-fly heuristic adaptation. We extend their approach to the pseudospectral approximation of functions. Adaptivity is expected to yield substantial efficiency gains in high dimensions—particularly for functions with anisotropic dependence on input parameters and functions whose inputs might not be strongly coupled at high order.

Previous attempts to extend pseudospectral methods to multivariate polynomial approximation with sparse model evaluations employed ad hoc approaches that are not always accurate. A common procedure has been to use sparse quadrature, or even dimension-adaptive sparse quadrature, to evaluate polynomial coefficients directly [39, 21]. This leads to at least two difficulties. First, the truncation of the polynomial expansion must be specified independently of the quadrature grid, yet it is unclear how to do this, particularly for anisotropic and generalized sparse grids. Second, unless one uses excessively high-order quadrature, significant aliasing errors may result. Constantine et al. [6] provided the first clear demonstration of these aliasing errors and proposed a Smolyak algorithm that does not share them. That work also demonstrated a link between Smolyak pseudospectral approximation and an extension to Lagrange interpolation called sparse interpolation, which uses function evaluations on a sparse grid and has well characterized convergence properties [24, 2].

The first half of this work performs a theoretical analysis, placing the solution from [6] in the broader context of Smolyak constructions, and explaining the origin of the observed aliasing errors for general (e.g., anisotropic) choices of sparse grid and quadrature rule. We do so by using the notion of polynomial exactness, without appealing to interpolation properties of particular quadrature rules. We establish conditions under which tensorized approximation operators are exact for particular polynomial inputs, then apply this analysis to the specific cases of quadrature and pseudospectral approximation; these cases are closely related and facilitate comparisons between Smolyak pseudospectral algorithms and direct quadrature. Section LABEL:sec:approximations develops computable one-dimensional and tensorized approximations for these settings. Section LABEL:sec:smolyak describes general Smolyak algorithms and their properties, yielding our principal theorem about the polynomial exactness of Smolyak approximations, and then applies these results to quadrature and pseudospectral approximation. Section LABEL:sec:comparison compares the Smolyak approach to conventional direct quadrature. Our error analysis of direct quadrature shows why the approach goes wrong and allows us to draw an important conclusion: in almost all cases, direct quadrature is not an appropriate method for constructing polynomial expansions and should be superseded by Smolyak pseudospectral methods.

These results provide a rigorous foundation for adaptivity, which is the second focus of this paper. Adaptivity makes it possible to harness the full flexibility of Smolyak algorithms in practical settings. Section LABEL:sec:adaptive introduces a fully adaptive algorithm for Smolyak pseudospectral approximation, which uses a single tolerance parameter to drive iterative refinement of both the polynomial approximation space and the corresponding collection of model evaluation points. As the adaptive method is largely heuristic, Section LABEL:sec:experiments demonstrates the benefits of this approach with numerical examples.

2 Full tensor approximations

Tensorization is a common approach for lifting one-dimensional operators to higher dimensions. Not only are tensor products computationally convenient, but they provide much useful structure for analysis. In this section, we develop some essential background for computable tensor approximations, then apply it to problems of (i) approximating integrals with numerical quadrature; and (ii) approximating projection onto polynomial spaces with pseudospectral methods. In particular, we are interested in analyzing the errors associated with these approximations and in establishing conditions under which the approximations are exact.

2.1 General setting

Consider a collection of one-dimensional linear operators , where indexes the operators used in different dimensions. In this work, will be either an integral operator or an orthogonal projector onto some polynomial space. We can extend a collection of these operators into higher dimensions by constructing the tensor product operator


The one-dimensional operators need not be identical; the properties of the resulting tensor operator are constructed independently from each dimension. The bold parenthetical superscript refers to the tensor operator instead of the constituent one-dimensional operators.

As the true operators are not available computationally, we work with a convergent sequence of computable approximations, , such that


in some appropriate norm. Taking the tensor product of these approximations provides an approximation to the full tensor operator, , where the level of the approximation may be individually selected in each dimension, so the tensor approximation is identified by a multi-index . Typically, and in the cases of interest in this work, the tensor approximation will converge in the same sense as the one-dimensional approximation as all components of .

An important property of approximation algorithms is whether they are exact for some inputs; characterizing this set of inputs allows us to make useful statements at finite order.

Definition 2.1 (Exact Sets)

For an operator and a corresponding approximation , define the exact set as and the half exact set .

The half exact set will help connect the exactness of a quadrature rule to that of the closely related pseudospectral operators. This notation is useful in proving the following lemma, which relates the exact sets of one-dimensional approximations and tensor approximations.

Lemma 2.2

If a tensor approximation is constructed from one-dimensional approximations with known exact sets, then


Proof. It is sufficient to show that the approximation is exact for an arbitrary monomial input where , because we may extend to sums by linearity:

The first step uses the tensor product structure of the operator and the second uses the definition of exact sets.     

2.2 Multi-indices

Before continuing, we must make a short diversion to multi-indices, which provide helpful notation when dealing with tensor problems. A multi-index is a vector . An important notion for multi-indices is that of a neighborhood.

Definition 2.3 (Neighborhoods of multi-indices)

A forward neighborhood of a multi-index is the multi-index set , where are the canonical unit vectors. The backward neighborhood of a multi-index is the multi-index set .

Smolyak algorithms rely on multi-index sets that are admissible.

Definition 2.4 (Admissible multi-indices and multi-index sets)

A multi-index is admissible with respect to a multi-index set if . A multi-index set is admissible if every is admissible with respect to .

Two common admissible multi-index sets with simple geometric structure are total order multi-index sets and full tensor multi-index sets. One often encounters total order sets in the sparse grids literature and full tensor sets when dealing with tensor grids of points. The total order multi-index set comprises those multi-indices that lie within a -dimensional simplex of side length :


The full tensor multi-index set is the complete grid of indices bounded term-wise by a multi-index :


2.3 Integrals and quadrature

Let be an open or closed interval of the real line . Then we define the weighted integral operator in one dimension as follows:


where is some real-valued function and is an integrable weight function. We may extend to higher dimensions by forming the tensor product integral , which uses separable weight functions and Cartesian product domains.

Numerical quadrature approximates the action of an integral operator with a weighted sum of point evaluations. For some family of quadrature rules, we write the “level ” quadrature rule, comprised of points, as


We call the growth rate of the quadrature rule, and its form depends on the quadrature family; some rules only exist for certain numbers of points and others may be tailored, for example, to produce linear or exponential growth in the number of quadrature points with respect to the level.

Many quadrature families are exact if is a polynomial of a degree or less, which allows us to specify a well-structured portion of the exact set for these quadrature rules:


where is the space of polynomials of degree or less. It is intuitive and useful to draw the exact set as in Figure LABEL:fig:Quadrature1D. For this work, we rely on quadrature rules that exhibit polynomial accuracy of increasing order, which is sufficient to demonstrate convergence for functions in [4].

Figure 2.1: Consider, one-dimensional Gaussian quadrature rule with three points, , which is exact for fifth degree polynomials. This diagram depicts the exact set, , and half exact set, , of this quadrature rule.

Tensor product quadrature rules are straightforward approximations of tensor product integrals that inherit convergence properties from the one-dimensional case. The exact set of a tensor product quadrature rule includes the tensor product of the constituent approximations’ exact sets, as guaranteed by Lemma LABEL:thm:tensorAccuracy and depicted in Figure LABEL:fig:SingleTensorQuadrature.

Figure 2.2: Consider a two dimensional quadrature rule constructed from three point Gaussian quadrature rules, . This diagram depicts the exact set, , and half exact set, .

2.4 Polynomial projection

A polynomial chaos expansion approximates a function with a weighted sum of orthonormal polynomials [38, 43]. Let be the separable Hilbert space of square-integrable functions , with inner product defined by the weighted integral , and normalized so that it may represent a probability density. Let be the space of univariate polynomials of degree or less. Now let be an orthogonal projector onto this subspace, written in terms of polynomials orthonormal with respect to the inner product of :


The polynomial space is of course the range of the projection operator. These polynomials are dense in , so the polynomial approximation of any converges in the sense as [4, 43]. If , the coefficients must satisfy .

Projections with finite degree omit terms of the infinite series, thus incurring truncation error. We can write this error as


Hence, we may reduce the truncation error to any desired level by increasing , removing terms from the sum in (LABEL:eq:truncation) [4, 17].

The -dimensional version of this problem requires approximating functions in the Hilbert space via a tensor product basis of the univariate polynomials defined above:


where . The multi-index tailors the range of the projection to include a rectangular subset of polynomials.

As in the one-dimensional case, truncation induces error equal to the sum of the squares of the omitted coefficients, which we may similarly reduce to zero as , . The multivariate polynomial expansion also converges in an sense for any [4]

2.5 Aliasing errors in pseudospectral approximation

The inner products defining the expansion coefficients above are not directly computable. Pseudospectral approximation provides a practical non-intrusive algorithm by approximating these inner products with quadrature. Define the pseudospectral approximation in one dimension as


where is the polynomial truncation at level , to be specified shortly [4, 17]. Pseudospectral approximations are constructed around a level quadrature rule, and are designed to include as many terms in the sum as possible while maintaining accuracy. Assuming that , we can compute the error between the pseudospectral approximation and an exact projection onto the same polynomial space:


This quantity is the aliasing error [4, 17]. The error is non-zero because quadrature in general only approximates integrals; hence each is an approximation of . The pseudospectral operator also incurs truncation error, as before, which is orthogonal to the aliasing error. We can expand each approximate coefficient as


The first step substitutes in the polynomial expansion of , which we assume is convergent, and rearranges using linearity. The second step partitions the sum around the truncation of the pseudospectral expansion. Although the basis functions are orthonormal, , we cannot assume in general that the approximation . Now substitute (LABEL:e:coeff) back into the aliasing error expression:


This form reveals the intimate link between the accuracy of pseudospectral approximations and the polynomial accuracy of quadrature rules. All aliasing is attributed to the inability of the quadrature rule to determine the orthonormality of the basis polynomials, causing one coefficient to corrupt another. The contribution of the first two parenthetical terms on the right of (LABEL:eqn:errorDecomp1D) is called internal aliasing, while the third term is called external aliasing. Internal aliasing is due to inaccuracies in when both and are included in the expansion, while external aliasing occurs when only one of these polynomials is included in the expansion. For many practical quadrature rules (and for all those used in this work), if and , and hence the discrete inner product is not zero, then is [34]. As a result, the magnitude of an aliasing error typically corresponds to the magnitude of the aliased coefficients.

In principle, both types of aliasing error are driven to zero by sufficiently powerful quadrature, but we are left to select for a particular quadrature level . External aliasing is in the magnitude of the truncation error, and thus it is driven to zero as long as increases with . Internal aliasing could be with respect to the function of interest, meaning that the procedure neither converges nor provides a useful approximation. Therefore, the obvious option is to include as many terms as possible while setting the internal aliasing to zero.

For quadrature rules with polynomial exactness, we may accomplish this by setting . This ensures that the internal aliasing of is zero, because . Equivalently, a pseudospectral operator using quadrature has a range corresponding to the half exact set . Alternatively, we may justify this choice by noting that it makes the pseudospectral approximation exact on its range, .

Given this choice of , we wish to show that the pseudospectral approximation converges to the true function, where the magnitude of the error is as follows:


The two terms on right hand side comprise the external aliasing and the truncation error, respectively. We already know that the truncation error goes to zero as . The external aliasing also vanishes for functions , as the truncated portion of likewise decreases [34]. In the case of Gaussian quadrature rules, a link to interpolation provides precise rates for the convergence of the pseudospectral operator based on the regularity of [4].

As with quadrature algorithms, our analysis of pseudospectral approximation in one dimension is directly extensible to multiple dimensions via full tensor products. We may thus conclude that converges to the projection onto the tensor product polynomial space in the same sense. The exact set follows Lemma LABEL:thm:tensorAccuracy, and hence the tensor product approximation inherits zero internal aliasing if suitable one-dimensional operators are used.

3 Smolyak algorithms

Thus far, we have developed polynomial approximations of multivariate functions by taking tensor products of one-dimensional pseudospectral operators. Smolyak algorithms avoid the exponential cost of full tensor products when the input dimensions are not fully coupled, allowing the use of a telescoping sum to blend different lower-order full tensor approximations.

Example 3.1

Suppose that . To construct a polynomial expansion with both the and terms, a full tensor pseudospectral algorithm would estimate all the polynomial terms up to , because tensor algorithms fully couple the dimensions. This mixed term is costly, requiring, in this case, an point grid for Gaussian quadratures. The individual terms can be had much more cheaply, using , , and grids, respectively. Smolyak algorithms help realize such savings in practice.

This section reviews the construction of Smolyak algorithms and presents a new theorem about the exactness of Smolyak algorithms built around arbitrary admissible index sets. We apply these results to quadrature and pseudospectral approximation, allowing a precise characterization of their errors.

3.1 General Smolyak algorithms

As in Section LABEL:sec:approximations, assume that we have for every dimension a convergent sequence of approximations. Let denote the collection of these sequences over all the dimensions. Define the difference operators


For any , we may write the exact or “true” operator as the telescoping series


Now we may write the tensor product of the exact operators as the tensor product of the telescoping sums, and interchange the product and sum:


Smolyak’s idea is to approximate the tensor product operator with truncations of this sum [31]:


We refer to the multi-index set as the Smolyak multi-index set, and it must be admissible for the sum to telescope correctly. Smolyak specifically suggested truncating with a total order multi-index set, which is the most widely studied choice. However, we can compute the approximation with any admissible multi-index set. Although the expression above is especially clean, it is not the most useful form for computation. We can reorganize the terms of (LABEL:eq:generalSmolyakDiff) to construct a weighted sum of the tensor operators:


where are integer Smolyak coefficients computed from the combinatorics of the difference formulation. One can compute the coefficients through a simple iteration over the index set and use (LABEL:eq:generalSmolyakDiff) to determine which full tensor rules are incremented or decremented. In general, these coefficients are non-zero near the leading surface of the Smolyak multi-index set, reflecting the mixing of the most accurate constituent full tensor approximations.

If each sequence of one-dimensional operators converges, then the Smolyak approximation converges to the tensor product of exact operators as . For the isotropic simplex index set, some precise rates of convergence are known with respect to the side length of the simplex [36, 37, 35, 29, 30]. Although general admissible Smolyak multi-index sets are difficult to study theoretically, they allow detailed customization to the anisotropy of a particular function.

3.2 Exactness of Smolyak algorithms

In the one-dimensional and full tensor settings, we have characterized approximation algorithms through their exact sets—those inputs for which the algorithm is precise. This section shows that if the constituent one-dimensional approximations have nested exact sets, Smolyak algorithms are the ideal blending of different full tensor approximations from the perspective of exact sets; that is, the exact set of the Smolyak algorithm contains the union of the exact sets of the component full tensor approximations. This result will facilitate subsequent analysis of sparse quadrature and pseudospectral approximation algorithms. This theorem and our proof closely follow the framework provided by Novak and Ritter [25, 26, 2], but include a generalization to arbitrary Smolyak multi-index sets.

Theorem 3.2

Let be a Smolyak algorithm composed of linear operators with nested exact sets, i.e., with implying that for , where is admissible. Then the exact set of contains


Proof. We begin by introducing notation to incrementally build a multi-index set dimension by dimension. For a multi-index set of dimension , let the restriction of the multi-indices to the first dimensions be . Furthermore, define subsets of based on the th element of the multi-indices, . These sets are nested, , because is admissible. Also let denote the maximum value of the th component of the multi-indices in the set .

Using this notation, one can construct inductively,


It is sufficient to prove that the Smolyak operator is exact for an arbitrary with tensor structure, . Suppose there exists a such that . We will show that if is an admissible multi-index set containing , then is exact on . We do so by induction on the dimension of the Smolyak operator and the function.

First, consider the case. , where . Hence

For the induction step, we construct the -dimensional Smolyak operator in terms of the -dimensional operator:


This sum is over increasing levels of accuracy in the dimension. We know the level required for the approximate operator to be exact in this dimension; this may be expressed as


Therefore the sum (LABEL:e:fullsum) can be truncated at the term, as the differences of higher terms are zero when applied to :


Naturally, . By nestedness, is also contained in for . The induction hypothesis then guarantees


Applying the -dimensional Smolyak operator to the truncated version of yields


Since each of the -dimensional Smolyak algorithms is exact, by the induction hypothesis, we replace them with the true operators and rearrange by linearity to obtain


The approximation in the dimension is exactly of the level needed to be exact on the th component of . Then (LABEL:e:lastdim) becomes


Thus the Smolyak operator is precise for , and the claim is proven.     

3.3 Smolyak quadrature

We recall the most familiar use of Smolyak algorithms, sparse quadrature. Consider a family of one-dimensional quadrature rules in each dimension ; denote these rules by . The resulting Smolyak algorithm is written as:


This approximation inherits its convergence from the one-dimensional operators. The set of functions that are exactly integrated by a Smolyak quadrature algorithm is described as a corollary of Theorem LABEL:thm:SmolyakAccuracy.

Corollary 3.3

For a sparse quadrature rule satisfying the hypotheses of Theorem LABEL:thm:SmolyakAccuracy,


Quadrature rules with polynomial accuracy do have nested exact sets, as required by the theorem. An example of Smolyak quadrature exact sets is shown in Figure LABEL:fig:smolyakQuadratureAccuracies.

(a) The exact set for a level-four Smolyak quadrature in two dimensions, based on linear growth Gaussian quadrature rules.
(b) The exact set for a level-three Smolyak quadrature in two dimensions, based on exponential growth Gaussian quadrature rules.
Figure 3.3: The exact set diagram for two Smolyak quadrature rules, and the corresponding basis for a Smolyak pseudospectral approximation. is shown in a solid line, is the dashed line. The staircase appearance results from the superposition of rectangular full tensor exact sets.

3.4 Smolyak pseudospectral approximation

Applying Smolyak’s algorithm to pseudospectral approximation operators yields a sparse algorithm that converges under similar conditions as the one-dimensional operators from which it is constructed. This algorithm is written as


The Smolyak algorithm is therefore a sum of different full tensor pseudospectral approximations, where each approximation is built around the polynomial accuracy of a single full tensor quadrature rule. It is not naturally expressed as a set of formulas for the polynomial coefficients, because different approximations include different polynomials. The term is included in the Smolyak approximation if and only if . Here, is the full tensor quadrature rule used by the full tensor pseudospectral approximation . As in the full tensor case, the half exact set of a Smolyak quadrature rule defines the range of the Smolyak pseudospectral approximation.

Once again, the Smolyak construction guarantees that the convergence of this approximation is inherited from its constituent one-dimensional approximations. Our choices for the pseudospectral operators ensure nestedness of the constituent exact sets, so we may use Theorem LABEL:thm:SmolyakAccuracy to ensure that Smolyak pseudospectral algorithms are exact on their range.

Corollary 3.4

If the constituent one-dimensional pseudospectral rules have no internal aliasing and satisfy the conditions of Theorem LABEL:thm:SmolyakAccuracy, then the resulting Smolyak pseudospectral algorithm has no internal aliasing.

We additionally provide a theorem that characterizes the external aliasing properties of Smolyak pseudospectral approximation, which the next section will contrast with direct quadrature.

Theorem 3.5

Let be a polynomial term included in the expansion provided by the Smolyak algorithm , and let be a polynomial term not included in the expansion. There is no external aliasing of onto if any of the following conditions is satisfied: (a) there exists a dimension for which ; or (b) there exists a multi-index such that is included in the range of and , where is the quadrature rule used in .

Proof. If condition (a) is satisfied, then and are orthogonal in dimension , and hence that inner product is zero. Every quadrature rule that computes the coefficient corresponding to basis term is accurate for polynomials of at least order . Since , every rule that computes the coefficient can numerically resolve the orthogonality, and therefore there is no aliasing. If condition (b) is satisfied, then the result follows from the cancellations exploited by the Smolyak algorithm, as seen in the proof of Theorem LABEL:thm:SmolyakAccuracy.     

These two statements yield extremely useful properties. First, any Smolyak pseudospectral algorithm, regardless of the admissible Smolyak multi-index set used, has no internal aliasing; this feature is important in practice and not obviously true. Second, while there is external aliasing as expected, the algorithm uses basis orthogonality to limit which external coefficients can alias onto an included coefficient. The Smolyak pseudospectral algorithm is thus a practically “useful” approximation, in that one can tailor it to perform a desired amount of work while guaranteeing reliable approximations of the selected coefficients. Computing an accurate approximation of the function only requires including sufficient terms so that the truncation and external aliasing errors are small.

4 Comparing direct quadrature to Smolyak pseudospectral approximation

The current UQ literature often suggests a direct quadrature approach for constructing polynomial chaos expansions [40, 21, 7, 18]. In this section, we describe this approach and show that, in comparison to a true Smolyak algorithm, it is inaccurate or inefficient in almost all cases. Our comparisons will contrast the theoretical error performance of the algorithms and provide simple numerical examples that illustrate typical errors and why they arise.

4.1 Direct quadrature polynomial expansions

At first glance, direct quadrature is quite simple. First, choose a multi-index set to define a truncated polynomial expansion:


The index set is typically admissible, but need not be. Second, select any -dimensional quadrature rule , and estimate every coefficient as:


Unlike the Smolyak approach, we are left to choose and independently, giving the appearance of flexibility. In practice, this produces a more complex and far more subtle version of the truncation trade-off discussed in Section LABEL:sec:approximations. Below, we will be interested in selecting and to replicate the quadrature points and output range of the Smolyak approach, as it provides a benchmark for achievable performance.

Direct quadrature does not converge for every choice of and ; consider the trivial case where does not grow infinitely. It is possible that including far too many terms in relative to the polynomial exactness of could lead to a non-convergent algorithm. Although this behavior contrasts with the straightforward convergence properties of Smolyak algorithms, most reasonable choices for direct quadrature do converge, and hence this is not our primary argument against the approach.

Instead, our primary concern is aliasing in direct quadrature and how it reduces performance at finite order. Both internal and external aliasing are governed by the same basic rule, which is just a restatement of how we defined aliasing in Section LABEL:sec:aliasing.

Remark 4.1

For a multi-index set and a quadrature rule , the corresponding direct quadrature polynomial expansion has no aliasing between two polynomial terms if .

The next two sections compare the internal and external aliasing with both theory and simple numeric examples.

4.2 Internal aliasing in direct quadrature

As an extension of Remark LABEL:rem:dqAliasing, direct quadrature has no internal aliasing whenever every pair has no aliasing. We can immediately conclude that for any basis set , there is some quadrature rule sufficiently powerful to avoid internal aliasing errors. In practice, however, this rule may not be a desirable one.

Example 4.2

Assume that for some function with two inputs, we wish to include the polynomial basis terms and . By Remark LABEL:rem:dqAliasing, the product of these two terms must be in the exact set; hence, the quadrature must include at least a full tensor rule of accuracy . Although we have not asked for any coupling, direct quadrature must assume full coupling of the problem in order to avoid internal aliasing.

Therefore we reach the surprising conclusion that direct quadrature inserts significant coupling into the problem, whereas we selected a Smolyak quadrature rule in hopes of leveraging the absence of that very coupling—making the choice inconsistent. For most sparse quadrature rules, we cannot include as many polynomial terms as the Smolyak pseudospectral approach without incurring internal aliasing, because the quadrature is not powerful enough in the mixed dimensions.

Example 4.3

Let be the two-dimensional domain . Select a uniform weight function, which corresponds to a Legendre polynomial basis. Let . Use Gauss-Legendre quadrature and an exponential growth rule, such that . Select a sparse quadrature rule based on a total order multi-index set . Figure LABEL:fig:ExpDirectQuadrature shows the exact set of this Smolyak quadrature rule (solid line) along with its half-exact set (dashed line), which encompasses all the terms in the direct quadrature polynomial expansion.

Now consider the polynomial, which is in the half-exact set. The product of the and polynomial terms is , which is not within the exact set of the sparse rule. Hence, aliases onto because this quadrature rule has limited accuracy in the mixed dimensions.

Using both the Smolyak pseudospectral and direct quadrature methods, we numerically compute the polynomial expansion for this example. The resulting coefficients are shown in Figure LABEL:fig:ExpInternalAliasing. Even though the two methods use the same information and project onto the same basis, the Smolyak result has no internal aliasing while direct quadrature shows significant internal aliasing. Although both methods correctly compute the coefficient, direct quadrature shows aliasing on as predicted, and also on , , and . In this case, direct quadrature is unable to determine the order of the input function or even whether the input is function of or . Alternating terms are computed correctly because of the parity of the functions.

The errors observed in this simple example demonstrate why it is crucial to eliminate internal aliasing in the construction of one-dimensional pseudospectral approximations and to ensure that the full tensor and Smolyak algorithms inherit that property. More complex functions demonstrate the same type of error, except that the errors resulting from multiple source terms are superimposed. Examples of the latter are given by Constantine et al. [6].

Figure 4.1: The exact set and polynomials included in the direct quadrature construction from Example LABEL:ex:directAliasing.
(a) Smolyak Pseudospectral
(b) Direct Quadrature
Figure 4.4: Numerical results for Example LABEL:ex:directAliasing; each color square indicates the log of the coefficient magnitude for the basis function at that position. The circle identifies the correct non-zero coefficient.

There are some choices for which direct quadrature has no internal aliasing: full tensor quadrature rules and, notably, Smolyak quadrature constructed from one-dimensional Gaussian quadrature rules with , truncated according to an isotropic total-order multi-index set. However, many useful sparser or more tailored Smolyak quadrature rules, e.g., based on exponential growth quadrature rules or adaptive anisotropic Smolyak index sets, will incur internal aliasing if the basis selection matches the range of the Smolyak algorithm. This makes them a poor choice when a comparable Smolyak pseudospectral algorithm uses the same evaluation points and produces an approximation with the same polynomial terms, but is guaranteed by construction to have zero internal aliasing. Alternately, it is possible to select a sufficiently small polynomial basis to avoid internal aliasing, but this approach requires unnecessary conservatism that could easily be avoided with a Smolyak pseudospectral approximation.

4.3 External aliasing

The difference in external aliasing between direct quadrature and Smolyak pseudospectral approximation is much less severe. Both methods exhibit external aliasing from terms far outside the range of the approximation, as such errors are a necessary consequence of using finite order quadrature. Since the methods are constructed from similar constituent one-dimensional quadrature rules, aliasing is of similar magnitude when it occurs.

Comparing Theorem LABEL:thm:smolyakExternal, condition (b), and Remark LABEL:rem:dqAliasing, we observe that if the direct quadrature method has no external aliasing between two basis terms, the equivalent Smolyak pseudospectral algorithm will not either. Yet the two methods perform differently because of their behavior on separable functions. Condition (a) of Theorem LABEL:thm:smolyakExternal provides an additional condition under which external aliasing will not occur under a Smolyak pseudospectral algorithm, and thus it has strictly less external aliasing in general.

Example 4.4

If we repeat Example LABEL:ex:directAliasing but choose to be a polynomial outside the approximation space, , we obtain the results in Figure LABEL:fig:externalAliasing. Now every non-zero coefficient is the result of external aliasing. Direct quadrature correctly computes some terms because of either parity or the few cases where Remark LABEL:rem:dqAliasing is satisfied. However, the Smolyak approach has fewer errors because the terms not between and are governed by condition (a) of Theorem LABEL:thm:smolyakExternal, and hence have no external aliasing.

This example is representative of the general case. Direct quadrature always incurs at least as much external aliasing as the Smolyak approach, and the methods become equivalent if the external term causing aliasing is of very high order. Although both methods will always exhibit external aliasing onto coefficients of the approximation for non-polynomial inputs, the truncation can in principle be chosen to include all the important terms, so that the remaining external aliasing is acceptably small.

(a) Smolyak Pseudospectral
(b) Direct Quadrature
Figure 4.7: Numerical results for Example LABEL:ex:externalAliasing; each color square indicates the log of the coefficient magnitude for the basis function at its position. The circle indicates the correct non-zero coefficient. The Smolyak pseudospectral approach has fewer terms corrupted by external aliasing in this case.

4.4 Summary of comparison

Compared to the Smolyak pseudospectral approach, direct quadrature yields larger internal and external aliasing errors. Because of these aliasing errors, direct quadrature is essentially unable to make efficient use of most sparse quadrature rules. The Smolyak pseudospectral approach, on the other hand, is guaranteed never to have internal aliasing if the one-dimensional pseudospectral operators are chosen according to simple guidelines. We therefore recommend against using direct quadrature. The remainder of the paper will focus on extensions of the basic Smolyak pseudospectral approach.

5 Adaptive polynomial approximations

When constructing a polynomial approximation of a black-box computational model, there are two essential questions: first, which basis terms should be included in the expansion; and second, what are the coefficients of those basis terms? The Smolyak construction allows detailed control over the truncation of the polynomial expansion and the work required to compute it. Since we typically do not have a priori information about the functional relationship generated by a black-box model, we develop an adaptive approach to tailor the Smolyak approximation to this function, following the dimension-adaptive quadrature approaches of Gerstner & Griebel [13] and Hegland [16].

The Smolyak algorithm is well suited to an adaptive approach. The telescoping sum converges in the same sense as the constituent one-dimensional operators as the index set grows to include , so we can simply add more terms to improve the approximation until we are satisfied. We separate our adaptive approach into two components: local improvements to the Smolyak multi-index set and a global stopping criterion.

5.1 Dimension adaptivity

Dimension adaptivity is responsible for identifying multi-indices to add to a Smolyak multi-index set in order to improve its accuracy. A standard refinement approach is simply to grow an isotropic simplex of side length . [13] and [16] instead suggest a series of greedy refinements that customize the Smolyak algorithm to a particular problem.

The refinement used in [13] is to select a multi-index and to add the forward neighbors of that are admissible. The multi-index is selected via an error indicator . We follow [13] and assume that whenever contributes strongly to the result of the algorithm, it represents a subspace that likely needs further refinement.

Let be a multi-index such that , where and are admissible multi-index sets. The triangle inequality (for some appropriate norm, see Section LABEL:s:adaptivecomments) bounds the change in the Smolyak approximation produced by adding to , yielding a useful error indicator:


Conveniently, this error indicator does not change as evolves, so we need only compute it once. At each adaptation step, we find the that maximizes and that has at least one admissible forward neighbor. Then we add those forward neighbors.

5.2 Termination criterion

Now that we have a strategy to locally improve a multi-index set, it is useful to have a global estimate of the error of the approximation, . We cannot expect to compute the exact error, but even a rough estimate is useful. We follow Gerstner & Griebel’s choice of global error indicator


where the sum is taken over all the multi-indices that are eligible for local adaptation at any particular step (i.e., that have admissible forward neighbors) [13]. The algorithm may be terminated when a particular threshold of the global indicator is reached, or when it falls by a specified amount.

5.3 Error indicators and work-considering algorithms

Thus far we have presented the adaptation strategy without reference to the problem of polynomial approximation. In this specific context, we use the norm in (LABEL:eq:localError), because it corresponds to the convergence properties of pseudospectral approximation and thus seems an appropriate target for greedy refinements. This choice is a heuristic to accelerate performance—albeit one that is simple and natural, and has enjoyed success in numerical experiments (see Section LABEL:sec:experiments). Moreover, the analysis of external aliasing in Theorem LABEL:thm:smolyakExternal suggests that, in the case of pseudospectral approximation, significant missing polynomial terms alias onto some of the included lower-order coefficients, giving the algorithm a useful indication of which direction to refine. This behavior helps reduce the need for smoothness in the coefficient pattern. Section LABEL:sec:quadRules provides a small fix that further helps with even or odd functions.

One is not required to use this norm to define , however, and it is possible that other choices could serve as better heuristics for some problems. Unfortunately, making definitive statements about the properties or general utility of these heuristic refinement schemes is challenging. The approach described above is intended to be broadly useful, but specific applications may require experimentation to find better choices.

Beyond the choice of norm, a commonly considered modification to is to incorporate a notion of the computational effort required to refine . Define as the amount of work to refine the admissible forward neighbors of , e.g., the number of new function evaluation points. [13] discusses an error indicator that provides a parameterized sliding scale between selecting the term with highest and the lowest :


Here is the tuning parameter, and and are the indicator and cost of the first term. Putting considers only the standard error indicator and considers only the cost. A different indicator with a similar intent is


where describes a conversion between error and work. Both of these methods will sometimes select terms of low cost even if they do not appear to provide immediate benefit to the approximation. Yet we find both methods to be challenging to use in practice, because of the difficulty of selecting the tuning parameter. One can remove this particular tuning requirement by taking a ratio:


This indicator looks for “efficient” terms to refine—ones that are expected to yield greater error reduction at less cost—rather than simply the highest-error directions. We performed some numerical experiments with these methods, but none of the examples demonstrated significant improvement. Furthermore, poor choices of tuning parameters can be harmful because they can essentially make the algorithm revert to a non-adaptive form. We do not give detailed results here for those experiments because they are not particularly conclusive; some types of coefficient patterns may benefit from work-considering approaches, but this remains an open problem.

On a similar note, using as a termination criteria is also a heuristic. As our experiments in Section LABEL:sec:experiments will show, for most smooth functions is an excellent estimate of the approximation accuracy. In other cases, the indicator can be quite poor; hence one should not rely on it exclusively. In practice, we typically terminate the algorithm based on a combination of elapsed wall clock time, the global error indicator, and an error estimate computed from limited ad hoc sampling.

6 Numerical experiments

Our numerical experiments focus on evaluating the performance of different quadrature rules embedded within the Smolyak pseudospectral scheme, and on evaluating performance of the adaptive Smolyak approximation strategy. Aside from the numerical examples of Section LABEL:sec:comparison, we do not investigate the performance of direct quadrature any further. Given our theoretical analysis of aliasing errors and the numerical demonstrations in [6], one can conclude without further demonstration that destructive internal aliasing indeed appears in practice.

This section begins by discussing practical considerations in the selection of quadrature rules. Then we evaluate convergence of Smolyak pseudospectral approximation schemes (non-adaptive and adaptive) on the Genz test functions. Next, we approximate a larger chemical kinetic system, illustrating the efficiency and accuracy of the adaptive method. Finally, we evaluate the quality of the global error indicator on all of these examples.

6.1 Selection of quadrature rules

Thus far we have sidestepped practical questions about which quadrature rules exist or are most efficient. Our analysis has relied only on polynomial accuracy of quadrature rules; all quadrature rules with a given polynomial accuracy allow the same truncation of a pseudospectral approximation. In practice, however, we care about the cumulative cost of the adaptive algorithm, which must step through successive levels of refinement.

Integration over a bounded interval with uniform weighting offers the widest variety of quadrature choices, and thus allows a thorough comparison. Table LABEL:tab:quadCost summarizes the costs of several common quadrature schemes. First, we see that linear-growth Gaussian quadrature is asymptotically much less efficient than exponential-growth in reaching any particular degree of exactness. However, for rules with fewer than about ten points, this difference is not yet significant. Second, Clenshaw-Curtis shows efficiency equivalent to exponential-growth Gaussian: both use points to reach th order polynomial exactness [5]. However, their performance with respect to external aliasing differs: Clenshaw-Curtis slowly loses accuracy if the integrand is of order greater than , while Gaussian quadrature gives error even on -order functions [34]. This may make Clenshaw-Curtis Smolyak pseudospectral estimates more efficient. Finally, we consider Gauss-Patterson quadrature, which is nested and has significantly higher polynomial exactness—for a given cumulative cost—than the other types [27]. Computing the quadrature points and weights in finite precision (even extended-precision) arithmetic has practically limited Gauss-Patterson rules to 255 points, but we recommend them whenever this is sufficient.

Lin. G Exp. G C-C G-P
1 1 1 1 1 1 1 1 1 1 1
2 2 3 3 2 3 3 3 3 3 5
3 3 5 6 4 7 7 5 5 7 10
4 4 7 10 8 15 15 9 9 15 22
5 5 9 15 16 31 31 17 17 31 46
6 6 11 21 32 63 63 31 31 63 94
Table 6.1: The cost of four quadrature strategies as their order increases: linear growth Gauss-Legendre quadrature (Lin. G), exponential growth Gauss-Legendre quadrature (Exp. G), Clenshaw-Curtis quadrature (C-C), and Gauss-Patterson quadrature (G-P). We list the number of points used to compute the given rule (p), the polynomial exactness (a), and the total number of points used so far (t). For nested rule, (p) = (t), so the total column is omitted.

For most other weights and intervals, there are fewer choices that provide polynomial exactness, so exponential-growth Gaussian quadrature is our default choice. In the specific case of Gaussian weight, Genz has provided a family of Kronrod extensions, similar to Gauss-Patterson quadrature, which may be a useful option [12].

If a linear growth rule is chosen and the domain is symmetric, we suggest that each new level include at least two points, so that the corresponding basis grows by at least one even and one odd basis function. This removes the possibility for unexpected effects on the adaptive strategy if the target function is actually even or odd.

6.2 Basic convergence: Genz functions

The Genz family [10, 11] comprises six parameterized functions, defined from . They are commonly used to investigate the accuracy of quadrature rules and interpolation schemes [2, 20]. The purpose of this example is to show that different Smolyak pseudospectral strategies behave roughly as expected, as evidenced by decreasing approximation errors as more function evaluations are employed. These functions are as follows:

Our first test uses five isotropic and non-adaptive pseudospectral approximation strategies. The initial strategy is the isotropic full tensor pseudospectral algorithm, based on Gauss-Legendre quadrature, with order growing exponentially with level. The other four strategies are total-order expansions of increasing order based on the following quadrature rules: linear growth Gauss-Legendre, exponential growth Gauss-Legendre, Clenshaw-Curtis, and Gauss-Patterson. All the rules were selected so that the final rule would have around points.

We consider 30 random realizations of each Genz function in dimensions; random parameters for the Genz functions are drawn uniformly from , then normalized so that and , where indexes the Genz function type and the constants are as chosen in [2, 20]. This experiment only uses the first four Genz functions, which are in , as pseudospectral methods have well known difficulties on functions with discontinuities or discontinuous derivatives [4]. Each estimate of approximation error is computed by Monte Carlo sampling with 10 samples. Figure LABEL:fig:GenzResults plots error at each stage, where each point represents the mean error over the 30 random functions.

Relatively simple conclusions can be drawn from this data. All the methods show fast convergence, indicating that the internal aliasing issues have indeed been resolved. In contrast, one would expect direct quadrature to suffer from large aliasing errors for the three super-linear growth rules. Otherwise, judging the efficiency of the different rules is not prudent, because differences in truncation and the structure of the test functions themselves obscure differences in efficiency. In deference to our adaptive strategy, we ultimately do not recommend this style of isotropic and function-independent truncation anyway.

To test our adaptive approach, Figure LABEL:fig:GenzScatter shows results from a similar experiment, now comparing the convergence of an adaptive Smolyak pseudospectral algorithm with that of a non-adaptive algorithm. To make the functions less isotropic, we introduce an exponential decay, replacing each with , where the are generated and normalized as above. For consistency, both algorithms are based on Gauss-Patterson quadrature. As we cannot synchronize the number of evaluations used by the adaptive algorithm for different functions, we plot individual errors for the 30 random functions instead of the mean error. This reveals the variability in difficulty of the functions, which was hidden in the previous plot. We conclude that the adaptive algorithm also converges as expected, with performance comparable to or better than the non-adaptive algorithm. Even though we have included some anisotropy, these functions include relatively high degrees of coupling; hence, in this case the non-adaptive strategy is a fairly suitable choice. For example, the “product peak” function shows little benefit from the adaptive strategy. Although omitted here for brevity, other quadrature rules produce similar results when comparing adaptive and non-adaptive algorithms.

Figure 6.1: Mean convergence of the non-adaptive isotropic total-order Smolyak pseudospectral algorithm with various quadrature rules, compared to the full tensor pseudospectral algorithm, on the Genz test functions.
Figure 6.2: convergence of the adaptive and non-adaptive Gauss-Patterson Smolyak pseudospectral algorithm. Individual results for 30 random instances of the Genz functions are shown.

6.3 Adaptivity: chemical kinetics

To further illustrate the benefits of an adaptive Smolyak approach, we build a surrogate for a realistic simulation of a combustion kinetics problem. Specifically, we consider the auto-ignition of a methane-air mixture given 14 uncertain rate parameters. Governing equations for this process are a set of stiff nonlinear ODEs expressing conservation of energy and of chemical species [19]. The uncertain rate parameters represent activation energies of reactions governing the conversion of methane to methyl, each endowed with a uniform distribution varying over of the nominal value. These parameters appear in Arrhenius expressions for the species production rates, with the reaction pathways and their nominal rate parameters given by the GRIMech 3.0 mechanism [15]. The output of interest is the logarithm of the ignition time, which is a functional of the trajectory of the ODE system, and is continuous over the selected parameter ranges. Simulations were performed with the help of the TChem software library [28], which provides convenient evaluations of thermodynamic properties and species production rates, along with Jacobians for implicit time integration.

Chemical kinetics are an excellent testbed for adaptive approximation because, by the nature of detailed kinetic systems, we expect strong coupling between some inputs and weak coupling between others, but we cannot predict these couplings a priori. We test the effectiveness of adaptive Smolyak pseudospectral methods based on the four quadrature rules discussed earlier. As our earlier analysis suggested that Gauss-Patterson quadrature should be most efficient, our reference solution is a non-adaptive Gauss-Patterson total-order Smolyak pseudospectral expansion. We ran the non-adaptive algorithm with a total order index set truncated at (which includes monomial basis terms up through ), using around 40000 point evaluations and taking over an hour to run. We tuned the four adaptive algorithms to terminate with approximately the same number of evaluations.

Figure LABEL:fig:combustionConvergence compares convergence of the five algorithms. The errors reported on the vertical axis are Monte Carlo estimates using points. Except for a small deviation at fewer than 200 model evaluations, all of the adaptive methods significantly outperform the non-adaptive method. The performance of the different quadrature rules is essentially as predicted in Section LABEL:sec:quadRules: Gauss-Patterson is the most efficient, exponential growth Gauss-Legendre and Clenshaw-Curtis are nearly equivalent, and linear growth Gauss-Legendre performs worse as the order of the polynomial approximation increases. Compared to the non-adaptive algorithm, adaptive Gauss-Patterson yields more than two orders of magnitude reduction in the error at the same number of model evaluations. Linear growth Gaussian quadrature is initially comparable to exponential growth Gaussian quadrature, because the asymptotic benefits of exponential growth do not appear while the algorithm is principally using very small one-dimensional quadrature rules. At the end of these experiments, a reasonable number of higher order quadrature rules are used and the difference becomes visible.

Figure 6.3: convergence of ignition delay in a 14-dimensional chemical kinetic system; comparing a non-adaptive isotropic total-order Gauss-Patterson-based Smolyak pseudospectral algorithm to the adaptive algorithm with various quadrature rules.
Figure 6.4: The plot depicts the difference between the number of coefficients of a particular magnitude and order in the final adaptive and non-adaptive Gauss-Patterson based expansions. The horizontal axis is the order of the term and the vertical axis specifies of the coefficient value. The color represents of the difference between the two methods, where positive values indicate more terms in the non-adaptive expansion. Hence, the dark blue at indicates that the non-adaptive expansion includes around 3,000 extra terms of magnitude and the dark red at indicates that the adaptive expansion includes about 1,000 extra terms of magnitude . Grey squares are the same for both expansions and white squares are not present in either.

We conclude by illustrating that the adaptive algorithm is effective because it successfully focuses its efforts on high-magnitude coefficients—that is, coefficients that make the most significant contributions to the function. Even though the non-adaptive expansion has around 37,000 terms and the final adaptive Gauss-Patterson expansion only has about 32,000 terms, the adaptive expansion exhibits much lower error because most of the additional terms in the non-adaptive expansion are nearly zero. By skipping many near-zero coefficients, the adaptive approach is able to locate and estimate a number of higher-order terms with large magnitudes. Figure LABEL:fig:combustionCoeffs depicts this pattern by plotting the difference between the numbers of included terms in the final adaptive Gauss-Patterson and non-adaptive expansions. The adaptive algorithm does not actually add any higher order monomials; neither uses one-dimensional basis terms of order higher than . Instead, the adaptive algorithm adds mixed terms of higher total order, thus capturing the coupling of certain variables in more detail than the non-adaptive algorithm. The figure shows that terms through 30th order are included in the adaptive expansion, all of which are products of non-constant polynomials in more than one dimension.

6.4 Performance of the global error indicator

To evaluate the termination criterion, we collected the global error indicator during runs of the adaptive algorithm for all of the test functions described above, including the slowly converging non-smooth Genz functions omitted before. The discontinuous Genz function does not include the exponential coefficient decay because the discontinuity already creates strong anisotropy. Results are shown for Gauss-Patterson quadrature. The relationship between the estimated error and the global error indicator is shown in Figure LABEL:fig:terminationPlot. For the smooth test functions, is actually an excellent indicator, as it is largely within an order of magnitude of the correct value and essentially linearly related to it. However, the non-smooth Genz functions illustrate the hazard of relying too heavily on this indicator: although the adaptive algorithm does decrease both the errors and the indicator, the relationship between the two appears far less direct.

Figure 6.5: Relationship between the termination criterion (LABEL:eq:globalError) and the estimated error for every function tested.

7 Conclusions

This paper gives a rigorous development of Smolyak pseudospectral algorithms, a practical approach for constructing polynomial chaos expansions from point evaluations of a function. A common alternative approach, direct quadrature, has previously been shown to suffer from large errors. We explain these errors as a consequence of internal aliasing and delineate the exact circumstances, derived from properties of the chosen polynomial basis and quadrature rules, under which internal aliasing will occur. Internal aliasing is a problem inherent to direct quadrature approaches, which use a single (sparse) quadrature rule to compute a set of spectral coefficients. These approaches fail because they substitute a numerical approximation for only a portion of the algorithm, i.e., the evaluation of integrals, without considering the impact of this approximation on the entire construction. For almost all sparse quadrature rules, internal aliasing errors may be overcome only through an inefficient use of function evaluations. In contrast, the Smolyak pseudospectral algorithm computes spectral coefficients by assembling tensor-product pseudospectral approximations in a coherent fashion that avoids internal aliasing by construction; moreover, it has smaller external aliasing errors. To establish these properties, we extend the known result that the exact set of a Smolyak pseudospectral approximation contains a union of the exact sets of all its constituent tensor-product approximation operators to the case of arbitrary admissible Smolyak multi-index sets. These results are applicable to any choice of quadrature rule and generalized sparse grid, and are verified through numerical demonstrations; hence, we suggest that the Smolyak pseudospectral algorithm is a superior approach in almost all contexts.

A key strength of Smolyak algorithms is that they are highly customizable through the choice of admissible multi-index sets. To this end, we describe a simple alteration to the adaptive sparse quadrature approaches of [13, 16], creating a corresponding method for adaptive pseudospectral approximation. Numerical experiments then evaluate the performance of different quadrature rules and of adaptive versus non-adaptive pseudospectral approximation. Tests of the adaptive method on a realistic chemical kinetics problem show multiple order-of-magnitude gains in accuracy over a non-adaptive approach. Although the adaptive strategy will not improve approximation performance for every function, we have little evidence that it is ever harmful and hence widely recommend its use.

While the adaptive approach illustrated here is deliberately simple, many extensions are possible. For instance, as described in Section LABEL:s:adaptivecomments, measures of computational cost may be added to the dimension refinement criterion. One could also use the gradient of the error indicator to identify optimal directions in the space of multi-indices along which to continue refinement, or to avoid adding all the forward neighbors of the multi-index selected for refinement. These and other developments will be pursued in future work.

A flexible open-source C++ code implementing the adaptive approximation method discussed in this paper is available at


The authors would like to thank Paul Constantine for helpful discussions and for sharing a preprint of his paper, which inspired this work. We would also like to thank Omar Knio and Justin Winokur for many helpful discussions, and Tom Coles for help with the chemical kinetics example. P. Conrad was supported during this work by a Department of Defense NDSEG Fellowship and an NSF graduate fellowship. P. Conrad and Y. Marzouk acknowledge additional support from the Scientific Discovery through Advanced Computing (SciDAC) program funded by the US Department of Energy, Office of Science, Advanced Scientific Computing Research under award number DE-SC0007099.


  • [1] I. Babuska, F. Nobile, and R. Tempone, A stochastic collocation method for elliptic partial differential equations with random input data, SIAM Journal on Numerical Analysis, 45 (2007), pp. 1005–1034.
  • [2] Volker Barthelmann, Erich Novak, and Klaus Ritter, High dimensional polynomial interpolation on sparse grids, Advances in Computational Mathematics, 12 (2000), pp. 273–288.
  • [3] John P Boyd, Chebyshev and Fourier Spectral Methods, Dover Publications, 2nd revise ed., 2001.
  • [4] Claudio Canuto, M. Yousuff Hussaini, Alfio Quarteroni, and Thomas A. Zhang Jr., Spectral Methods: Fundamentals in single domains, Springer Berlin Heidelberg, 2006.
  • [5] C.W. Clenshaw and A.R. Curtis, A method for numerical integration on an automatic computer, Numerische Mathematik, 2 (1960), pp. 197–205.
  • [6] Paul G Constantine, Michael S Eldred, and Eric T Phipps, Sparse Pseudospectral Approximation Method, Computer Methods in Applied Mechanics and Engineering, 229-232 (2012), pp. 1–12.
  • [7] M. S. Eldred and J. Burkardt, Comparison of non-intrusive polynomial chaos and stochastic collocation methods for uncertainty quantification, AIAA paper 2009–0976, (2009).
  • [8] Oliver G. Ernst, Antje Mugler, Hans-Jörg Starkloff, and Elisabeth Ullmann, On the convergence of generalized polynomial chaos expansions, ESAIM: Mathematical Modeling and Numerical Analysis, 46 (2012), pp. 317–339.
  • [9] W. Gander, Change of basis in polynomial interpolation, Numerical Linear Algebra with Applications, 12 (2005), pp. 769–778.
  • [10] Alan Genz, Testing Multidimensional Integration Routines, in Tools, Methods, and Languages for Scientific and Engineering Computation, B. Ford, J.C. Rault, and F. Thomasset, eds., North-Holland, 1984, pp. 81–94.
  • [11]  , A Package for Testing Multiple Integration Subroutines, in Numerical Integration: Recent Developments, Software and Applications, Patrick Keast and Graeme Fairweather, eds., D Reidel, 1987, pp. 337–340.
  • [12]  , Fully symmetric interpolatory rules for multiple integrals over infinite regions with Gaussian weight, Journal of Computational and Applied Mathematics, 71 (1996), pp. 299–309.
  • [13] T. Gerstner and M. Griebel, Dimension-Adaptive Tensor-Product Quadrature, Computing, 71 (2003), pp. 65–87.
  • [14] Roger Ghanem and Pol Spanos, Stochastic finite elements: a spectral approach, Springer-Verlag, 1991.
  • [15] Zhiwei Qin Gregory P. Smith, David M. Golden, Michael Frenklach, Nigel W. Moriarty, Boris Eiteneer, Mikhail Goldenberg, C. Thomas Bowman, Ronald K. Hanson, Soonho Song, William C. Gardiner, Jr., Vitali V. Lissianski, GRI-MECH 3.0.
  • [16] M Hegland, Adaptive sparse grids, ANZIAM J, 44 (2003), pp. 335–353.
  • [17] Jan S. Hesthaven, Sigal Gottlieb, and David Gottlieb, Spectral Methods for Time-Dependent Problems, Cambridge University Press, 2007.