Non-intrusive uncertainty quantification using reduced cubature rules

TL;DR

This paper introduces a novel method for generating nested quadrature and cubature rules with positive weights, enhancing uncertainty quantification by reducing nodes while maintaining accuracy and positivity, applicable in high-dimensional problems.

Contribution

The paper presents a reduction procedure for creating nested quadrature rules with positive weights, applicable to multi-dimensional sparse rules, improving upon existing methods like Smolyak rules.

Findings

Generated sparse cubature rules with positive weights.

Achieved competitive accuracy with fewer nodes.

Demonstrated effectiveness in fluid dynamics and mathematical test problems.

Abstract

For the purpose of uncertainty quantification with collocation, a method is proposed for generating families of one-dimensional nested quadrature rules with positive weights and symmetric nodes. This is achieved through a reduction procedure: we start with a high-degree quadrature rule with positive weights and remove nodes while preserving symmetry and positivity. This is shown to be always possible, by a lemma depending primarily on Carath\'eodory's theorem. The resulting one-dimensional rules can be used within a Smolyak procedure to produce sparse multi-dimensional rules, but weight positivity is lost then. As a remedy, the reduction procedure is directly applied to multi-dimensional tensor-product cubature rules. This allows to produce a family of sparse cubature rules with positive weights, competitive with Smolyak rules. Finally the positivity constraint is relaxed to allow more flexibility in the removal of nodes. This gives a second family of sparse cubature rules, in which iteratively as many nodes as possible are removed. The new quadrature and cubature rules are applied to test problems from mathematics and fluid dynamics. Their performance is compared with that of the tensor-product and standard Clenshaw-Curtis Smolyak cubature rule.

Tables

Table 1. Table 1 : Number of nodes of several cubature rules for several dimensions ( d 𝑑 d ) and several degrees ( K 𝐾 K ). N Positive subscript 𝑁 Positive N_{\textrm{Positive}} denotes the number of nodes of the symmetric reduced cubature rule with positive weights, N Negative subscript 𝑁 Negative N_{\textrm{Negative}} denotes the number of nodes of the negative symmetric reduced cubature rule, and N Smolyak subscript 𝑁 Smolyak N_{\textrm{Smolyak}} denotes the number of nodes of the smallest Smolyak sparse grid of at least degree K 𝐾 K .

$d$	$K$	$dim\mathbb{P}(K,d)$	$N_{\text{Smolyak}}$	$N_{\text{Positive}}$	$N_{\text{Negative}}$
5	5	252	61	113	43
5	7	792	241	544	384
5	9	2 002	805	1 313	325
5	11	4 368	2 473	4 096	2 016
5	13	8 568	7 245	6 005	1 607
7	5	792	113	689	99
7	7	3 432	589	1 797	325
7	9	11 440	2 471	19 717	901
7	11	31 824	9 101	28 479	2 863
7	13	77 520	30 907	158 709	28 479
10	5	3 003	221	13 461	201
10	7	19 448	1 581	20 533	1 361
10	9	92 378	8 810	1 368 449	3 705
10	11	352 716	41 445	8 284 617	12 489
10	13	1 144 066	172 055	26 598 325	38 353
15	5	15 504	481		451
15	9	1 307 504	40 001		30 861
15	13	37 442 160	1 472 697		362 063
20	5	53 130	841		801
20	9	10 015 005	120 401		98 881
25	5	142 506	1 301		1 251
25	9	52 451 256	286 001		244 101

Table 2. Table 2 : The test functions from Genz [ 8 ] . All functions are from a certain integrand family and depend on the parameters 𝐚 = ( a 1 , … , a N ) T 𝐚 superscript subscript 𝑎 1 … subscript 𝑎 𝑁 T \mathbf{a}={(a_{1},\dots,a_{N})}^{\mathrm{T}} and 𝐮 = ( u 1 , … , u N ) T 𝐮 superscript subscript 𝑢 1 … subscript 𝑢 𝑁 T \mathbf{u}={(u_{1},\dots,u_{N})}^{\mathrm{T}} . The parameter 𝐮 𝐮 \mathbf{u} is a parameter that does not affect the difficulty of the integral. The parameter 𝐚 𝐚 \mathbf{a} determines the degree to which the family attribute is present.

Integrand Family	Attribute
$f_1\left(𝐱\right)=\cos \left(2\pi u_1+{\sum }_{i=1}^n{a}_i{x}_i\right)$	Oscillatory
$f_2\left(𝐱\right)={\prod }_{i=1}^n{\left(a_i^{-2}+{\left(x_i-u_i\right)}^2\right)}^{-1}$	Product Peak
$f_3\left(𝐱\right)={\left(1+{\sum }_{i=1}^n{a}_i{x}_i\right)}^{-\left(n+1\right)}$	Corner Peak
$f_4\left(𝐱\right)=\exp \left(-{\sum }_{i=1}^n{a}_i^2{\left(x_i-u_i\right)}^2\right)$	Gaussian
$f_5\left(𝐱\right)=\exp \left(-{\sum }_{i=1}^n{a}_i\left|x_i-u_i\right|\right)$	$C_0$ function
	Discontinuous

Table 3. Table 3 : Uncertain parameters and their distribution as considered for the lid-driven cavity flow problem.

Parameter	Distribution
$u_{\text{lid}}$ (speed of the lid)	$\beta (3,3)$ with range $(0.5,1.5)$
$\nu$ (viscosity)	$\beta (4,4)$ with range $(0.0038,0.05)$

Table 4. Table 4 : Uncertain parameters and their distribution as considered for the aircraft aerodynamics test case.

Parameter	Distribution
Leading edge radius	$𝒩$ with 5% standard deviation
Maximum camber as percentage of the chord	$𝒩$ with 5% standard deviation
Distance of maximum camber from leading edge	$𝒩$ with 5% standard deviation
Angle of incidence	$\beta (4,4)$ with range ${2.31}^{\circ }\pm 5\%$
Side-slip angle	$\beta (4,4)$ with range $0^{\circ }\pm {0.5}^{\circ }$
Mach number	$\beta (4,4)$ with range $0.72\pm 5\%$
Free-stream pressure	$\beta (4,4)$ with range $\mathrm{101\hspace{0.17em}325}\text{N}/{\text{m}}^2\pm 5\%$

Table 5. Table 5 : The first four non-central moments determined either using the cubature rule directly on the results (without hat) or using a high-degree cubature rule on the least-squares estimation (with hat). Empty places are values smaller than 10 − 5 superscript 10 5 10^{-5} .

#	${𝐜}_{𝐥}$	${\widehat{𝐜}}_{𝐥}$	${𝐜}_{𝐝}$	${\widehat{𝐜}}_{𝐝}$	${𝐜}_{\text{sf}}$	${\widehat{𝐜}}_{\text{sf}}$
0	1.0000	1.0000	1.0000	1.0000	1.0000	1.0000
1	0.3640	0.3645	0.0226	0.0228	-0.0015	-0.0002
2	0.1326	0.1330	0.0005	0.0005
3	0.0483	0.0485
4	0.0176	0.0177

Equations

$$v\coloneqq v\bigl{(}s({\boldsymbol{\xi}})\bigr{)},$$

$$R(s;\mathbf{\xi })=0,$$

$$\mathbb{E}[u^l(\mathbf{\xi })]:={\int }_{\Xi }{u}^l(\mathbf{\xi })\mathrm{d}P(\mathbf{\xi }),\text{ for }l=1,2,\dots$$

$$\mathbb{E}[u^l(\mathbf{\xi })]\simeq \sum _{k=1}^N{u}^l({\mathbf{\xi }}_k)w_k,\text{ for }l=1,2,\dots ,$$

$$R(s({\mathbf{\xi }}_k);{\mathbf{\xi }}_k)=0$$

$$\sum _{k=1}^N{\xi }_k^j{w}_k={\int }_{\Xi }{\xi }^j\mathrm{d}P(\xi ),\text{ for all }j=0,\dots ,N-1.$$

$$\sum _{k=1}^N{m}_j({\mathbf{\xi }}_k)w_k={\int }_{\Xi }{m}_j(\mathbf{\xi })\mathrm{d}P(\mathbf{\xi }),\text{ for all }j=1,\dots ,N,$$

$${\mathcal{S}}_K=\sum _{\begin{array}{c}K-d+1\le \Vert \alpha {\Vert }_1\le K\\ \alpha \in {\mathbb{N}}^d\end{array}}{(-1)}^{K-\Vert \alpha {\Vert }_1}\left(\genfrac{}{}{0ex}{}{d-1}{K-\Vert \alpha {\Vert }_1}\right)\underset{k=1}{\overset{d}{⨂}}{\mathcal{Q}}_{N_{{\alpha }_k}},$$

$$N_k=\{\begin{array}{ll}1& \text{if k=1,}\\ 2^{k-1}+1& \text{otherwise.}\end{array}$$

$$\kappa :=\frac{{\sum }_{k=1}^N|w_k|}{{\sum }_{k=1}^N{w}_k}=\sum _{k=1}^N|w_k|=\mathcal{O}\left((\log N{)}^{d-1}\right).$$

$$\underset{V}{\underbrace{\left(\begin{array}{cccc}{\xi }_1^0& {\xi }_2^0& \dots & {\xi }_N^0\\ {\xi }_1^1& {\xi }_2^1& \dots & {\xi }_N^1\\ ⋮& ⋮& \ddots & ⋮\\ {\xi }_1^{N-2}& {\xi }_2^{N-2}& \dots & {\xi }_N^{N-2}\\ {\xi }_1^{N-1}& {\xi }_2^{N-1}& \dots & {\xi }_N^{N-1}\end{array}\right)}}\left(\begin{array}{c}{w}_1\\ w_2\\ ⋮\\ w_{N-1}\\ w_N\end{array}\right)=\left(\begin{array}{c}{\int }_{\Xi }{\xi }^0\mathrm{d}P(\xi )\\ {\int }_{\Xi }{\xi }^1\mathrm{d}P(\xi )\\ ⋮\\ {\int }_{\Xi }{\xi }^{N-2}\mathrm{d}P(\xi )\\ {\int }_{\Xi }{\xi }^{N-1}\mathrm{d}P(\xi )\end{array}\right),$$

$$\underset{V_{-1}}{\underbrace{\left(\begin{array}{cccc}{\xi }_1^0& {\xi }_2^0& \dots & {\xi }_N^0\\ {\xi }_1^1& {\xi }_2^1& \dots & {\xi }_N^1\\ ⋮& ⋮& \ddots & ⋮\\ {\xi }_1^{N-2}& {\xi }_2^{N-2}& \dots & {\xi }_N^{N-2}\end{array}\right)}}\left(\begin{array}{c}{w}_1\\ w_2\\ ⋮\\ w_{N-1}\\ w_N\end{array}\right)=\left(\begin{array}{c}{\int }_{\Xi }{\xi }^0\mathrm{d}P(\xi )\\ {\int }_{\Xi }{\xi }^1\mathrm{d}P(\xi )\\ ⋮\\ {\int }_{\Xi }{\xi }^{N-2}\mathrm{d}P(\xi )\end{array}\right),$$

$$\sum _{k=1}^{N+1}{c}_k{\mathbf{v}}_k=0.$$

$$\mathbf{v}$$

$$=\sum _{k=1}^{N+1}({\lambda }_k-\alpha c_k){\mathbf{v}}_k.$$

$$\alpha =\underset{k=1,\dots ,N+1}{\min }\left\{\frac{{\lambda }_k}{c_k}:c_k>0\right\}=:\frac{{\lambda }_{k_0}}{c_{k_0}}.$$

$$\mathbf{v}$$

$$\{{\xi }_1,{\xi }_2,\dots ,{\xi }_N\}=\{{\xi }_1,{\xi }_2,\dots ,{\xi }_{\frac{N}{2}},-{\xi }_{\frac{N}{2}},\dots ,-{\xi }_2,-{\xi }_1\}.$$

$$V_{-2}^{\prime }=\left(\begin{array}{cccc}{\xi }_1^0+{\xi }_N^0& {\xi }_2^0+{\xi }_{N-1}^0& \dots & {\xi }_{\frac{N}{2}}^0+{\xi }_{N-\frac{N}{2}+1}^0\\ {\xi }_1^1+{\xi }_N^1& {\xi }_2^1+{\xi }_{N-1}^1& \dots & {\xi }_{\frac{N}{2}}^1+{\xi }_{N-\frac{N}{2}+1}^1\\ {\xi }_1^2+{\xi }_N^2& {\xi }_2^2+{\xi }_{N-1}^2& \dots & {\xi }_{\frac{N}{2}}^2+{\xi }_{N-\frac{N}{2}+1}^2\\ ⋮& ⋮& \ddots & ⋮\\ {\xi }_1^{N-3}+{\xi }_N^{N-3}& {\xi }_2^{N-3}+{\xi }_{N-1}^{N-3}& \dots & {\xi }_{\frac{N}{2}}^{N-3}+{\xi }_{N-\frac{N}{2}+1}^{N-3}\end{array}\right).$$

$${\xi }_k^p+{\xi }_{N-k+1}^p={\xi }_k^p+(-1{)}^p{\xi }_k^p=0.$$

$$A_{-2}=\left(\begin{array}{cccc}2{\xi }_1^0& 2{\xi }_2^0& \dots & 2{\xi }_{\frac{N}{2}}^0\\ 2{\xi }_1^2& 2{\xi }_2^2& \dots & 2{\xi }_{\frac{N}{2}}^2\\ ⋮& ⋮& \ddots & ⋮\\ 2{\xi }_1^{N-4}& 2{\xi }_2^{N-4}& \dots & 2{\xi }_{\frac{N}{2}}^{N-4}\end{array}\right).$$

$$\{{\xi }_1,{\xi }_2,\dots ,{\xi }_N\}=\{{\xi }_1,{\xi }_2,\dots ,{\xi }_{\lfloor \frac{N}{2}\rfloor },0,-{\xi }_{\lfloor \frac{N}{2}\rfloor },\dots ,-{\xi }_1\}.$$

$$A_{-2}=\left(\begin{array}{ccccc}2{\xi }_1^0& 2{\xi }_2^0& \dots & 2{\xi }_{\lfloor \frac{N}{2}\rfloor }^0& 1\\ 2{\xi }_1^2& 2{\xi }_2^2& \dots & 2{\xi }_{\lfloor \frac{N}{2}\rfloor }^2& 0\\ ⋮& ⋮& \ddots & ⋮& ⋮\\ 2{\xi }_1^{N-3}& 2{\xi }_2^{N-3}& \dots & 2{\xi }_{\lfloor \frac{N}{2}\rfloor }^{N-3}& 0\end{array}\right).$$

$$\left(\genfrac{}{}{0ex}{}{\lfloor \frac{K}{2}\rfloor +d}{d}\right)<N_{\bar{Q}}.$$

$$N_Q=N_{\bar{Q}},$$

$$N_Q-\left(\genfrac{}{}{0ex}{}{\lfloor \frac{K-1}{2}\rfloor +d}{d}\right).$$

$$\left(\genfrac{}{}{0ex}{}{\lfloor \frac{K-1}{2}\rfloor +d}{d}\right).$$

$$1+\sum _{l=1}^K{p}_d(l)<N_{\bar{Q}},$$

$$1+\sum _{l=1}^{\lfloor K/2\rfloor }{p}_d(l)<N_{\bar{Q}},$$

$$p(x;a,b)\propto x^{a-1}(1-x{)}^{b-1}\text{ for }0\le x\le 1.$$

$$ϵ(N):=\Vert \bar{u}(N)-{\bar{u}}^{\ast }{\Vert }_2.$$

Full text

Keywords: Uncertainty Quantification, Numerical Integration, Cubature Rules

1 Introduction

The problem of non-intrusive uncertainty quantification (UQ) in expensive computational models is considered, for example computational fluid dynamics (CFD) models. Consider a model with $d$ uncertain parameters having specified distributions. The objective is to obtain statistics on the outputs of the model, while using the model only as a black box (i.e. non-intrusively). We wish to obtain accurate statistics with as few evaluations of the model as possible.

The canonical method is Monte Carlo (MC), with the well-known dimension independent convergence rate of $\mathcal{O}(1/\sqrt{N})$, where $N$ is the number of samples. For sufficiently low dimension $d$ this can be improved to \mathcal{O}\bigl{(}(\log N)^{d}/N\bigr{)} using Quasi Monte Carlo methods, see e.g. [32, 2]. For $d\lesssim 10$ this can be significantly further improved by using methods based on polynomial approximation of the model output in the parameter space. This case is studied in this paper. For a sufficiently smooth parametrized model, spectral convergence is obtained. Stochastic Collocation (SC) [14, 5, 31] is such a method which uses either tensor products or sparse grids to sample the parameter space. Quadrature weights on these grids allow the evaluation of statistics. We mention also the hybrid techniques of Witteveen et al. [27, 28], which use piecewise polynomial interpolation on random MC grids. Approaches based on compressed sensing (and therefore not dependent on quadrature rules) have also been studied [4, 1].

SC methods can be regarded as cubature rules targeted at moderate dimensional spaces. The conventional tensor-product cubature rule introduces a large number of nodes for moderate $d$. Sparse grid strategies are therefore required, e.g. a Smolyak sparse grid [24, 19]. These in turn require nested one-dimensional quadrature rules for the optimal result, but no general strategy exists to create these nested quadrature rules with positive weights for arbitrary distributions. Furthermore, a Smolyak procedure does not guarantee that the weights of the resulting multi-dimensional rule are positive (and therefore is not necessarily numerically stable), even if the underlying one-dimensional rule has positive weights.

Sparse grid techniques have been studied thoroughly by various authors. For example, Garcke, Gerstner, and Griebel (see e.g. [7, 9]) studied the generation of sparse grids, among others with improvements such as dimension dependent adaptivity. Narayan and Jakeman [15] studied the construction of quadrature rules as input for the Smolyak sparse grid. Nobile et al. [16] studied the effectiveness of sparse grids compared with MC methods. Anisotropic extensions (i.e. different quadrature rules in different dimensions) were also studied by Nobile et al. [17]. Pflüger [22] studied adaptive sparse grids, where locally the grid is refined if necessary, yielding a strategy to determine a sparse grid that depends on the specifics of the model.

The present paper has two major contributions. Firstly a method is introduced for constructing a nested family of one-dimensional quadrature rules with positive weights, from any single high-order rule (with positive weights). Thus given a quadrature rule for a specific probability distribution, a nested family can be constructed, suitable for use in a Smolyak procedure. If the original rule is symmetric, the symmetry of the entire family is guaranteed. Secondly, in the multi-dimensional case a closely related operation can be performed. Starting from a tensor-product rule, nodes can be removed successively while maintaining positivity of all weights and symmetry of the rule. The result is a new kind of sparse grid with only positive weights. Because the positivity restriction is quite limiting, also the case where negative weights are permitted is considered. This allows to remove more nodes at each step of the reduction procedure. All resulting rules are well suited for UQ, as top level quadrature rules can be chosen separately for each parameter, without any concern about nesting. The new rules are demonstrated on the Genz test functions, two CFD test cases, and compared to the tensor-product and Smolyak rules.

The study is set up as follows. First, in the section hereafter the UQ problem is formulated. In the next section some useful well-known methods are discussed. In Section 4 the reduced quadrature rule is introduced, which is extended to a multi-dimensional setting in Section 5. The introduced cubature rule is compared with conventional cubature rules in Section 6. Firstly a mathematical comparison is made using test functions. Secondly the cubature rules are applied in UQ for the standard lid-driven cavity flow problem computed through a Lattice Boltzmann method with two uncertain parameters. To show the effectiveness in high-dimensional problems, it is finally applied to a three-dimensional aircraft aerodynamics problem, computed through a finite-volume Euler-flow model, considering seven uncertain parameters.

2 Uncertainty Quantification

Consider a discrete computational problem for a quantity of interest

[TABLE]

where $s:\mathbb{R}^{d}\rightarrow\mathbb{R}^{n}$ is the state of some system, satisfying

[TABLE]

where $R$ is typically a discretization of a continuous PDE, including initial and boundary conditions, and where $n$ is the dimension of the discrete state. The quantity of interest $v:\mathbb{R}^{n}\rightarrow\mathbb{R}$ is a single quantity derived from the full state. The parameters ${\boldsymbol{\xi}}$ are $d$ random variables, that is ${\boldsymbol{\xi}}:\Omega\rightarrow\Xi$, which are assumed to be independent and square-integrable (i.e. having finite variance), with respect to the probability space $(\Omega,\mathcal{F},P)$ with $\Xi\subset\mathbb{R}^{d}$, $\Omega\subset\mathbb{R}^{d}$, $\mathcal{F}\subset 2^{\Omega}$, and $P$ the probability measure. Although infinite dimensional, random fields can be fit into this framework after the application of a truncated Karhunen-Loève expansion [10, 30].

The problem is now to determine the probability distribution and statistical moments of $u$, with u({\boldsymbol{\xi}})\coloneqq v\bigl{(}s({\boldsymbol{\xi}})\bigr{)}. The focus is on the latter, i.e. on determining

[TABLE]

The collocation approach is to approximate this integral using a weighted combination of a finite number of samples $\{\boldsymbol{\xi}_{k}\}_{k=1,\dots,N}\in\Xi$ as

[TABLE]

where $\{w_{k}\}_{k=1,\dots,N}\in\mathbb{R}$ are the weights. $u^{l}(\boldsymbol{\xi}_{k})=(u(\boldsymbol{\xi}_{k}))^{l}$ is determined by solving the (potentially expensive) deterministic discrete problem

[TABLE]

for $s({\boldsymbol{\xi}}_{k})$, and by evaluating $u({\boldsymbol{\xi}}_{k})$.

In the remainder of this paper the term quadrature rule is used in a one-dimensional setting (i.e. $d=1$) and the term cubature rule otherwise. All properties of cubature rules also apply to quadrature rules (but not vice versa).

3 Numerical integration - Terminology and basic principles

3.1 Quadrature and cubature rules

Let $\mathbb{P}(K,d)$ be all $d$-variate polynomials of degree equal to or less than $K$. The degree of a cubature rule is defined as the number $K$ such that all polynomials $p\in\mathbb{P}(K,d)$ are integrated exactly and at least one polynomial $p\in\mathbb{P}(K+1,d)$ exists that is not integrated exactly.

We consider a set of cubature rules to be nested if the nodes of a smaller cubature rule are also nodes of all larger cubature rules. If a cubature rule is nested, error estimates can be naturally constructed by comparing the approximation on two consecutive levels. In addition to nesting, it is desirable that rules are (i) symmetric, meaning the nodes and weights have the same symmetry as the underlying probability distribution, and (ii) positive, meaning all weights are positive. Symmetric quadrature rules naturally represent the underlying distribution and are necessary in the multi-dimensional case to reduce the number of nodes (which is done in the second part of this paper). Quadrature rules with positive weights are unconditionally numerically stable if they are evaluated and yield an integration operator with norm equal to 1. For example Gaussian quadrature rules are positive, irrespective of the underlying distribution, and symmetric if the distribution is symmetric [12]. However they are not nested. The nested Clenshaw–Curtis rule is usually applied with a uniform distribution, in which case weights are positive [3] – but this is not true if weights are constructed for an arbitrary distribution.

3.2 The generalized Vandermonde-matrix

If $N$ distinct one-dimensional quadrature nodes (denoted by $\{\xi_{k}\}_{k=1}^{N}\subset\mathbb{R}$) are specified, the weights of the quadrature rule can be determined such that it is a rule of degree $N-1$ by solving the following linear system:

[TABLE]

The matrix of this system is a Vandermonde-matrix, hereafter denoted by $V$ and defined by $V_{j,k}=\xi_{k}^{j}$. This system has a unique solution for distinct nodes. Hence the quadrature rule is unambiguously specified by the nodes only. For large $N$ the Vandermonde-matrix becomes ill-conditioned, such that for various quadrature rules more efficient algorithms exist to determine both the nodes and the weights, e.g. the algorithm of Golub and Welsch [12] can be used to determine a Gauss quadrature rule and Clenshaw–Curtis rules can be determined efficiently using a Fast Fourier transform [25].

Generalizing to a multi-dimensional setting, let $\{{\boldsymbol{\xi}}_{k}\}_{k=1}^{N}\subset\mathbb{R}^{d}$ be $N$ cubature nodes. Integration conditions result in the system

[TABLE]

where $m_{j}$ is the $j^{\text{th}}$ monomial under some ordering. We call the matrix $G_{j,k}=m_{j}(\boldsymbol{\xi}_{k})$ the generalized Vandermonde-matrix. As is well-known $G$ may be singular, but for tensor-product rules $G$ is non-singular, as it can be formed by the Kronecker product of the Vandermonde-matrix of the quadrature rules [13]. Just as in the one-dimensional case, in general this matrix can become ill-conditioned for high-dimensional polynomial spaces or large $N$.

3.3 Smolyak cubature rules

The Smolyak procedure [24] is a method of constructing “sparse” cubature rules from a family of (typically nested) quadrature rules indexed by level. Rather than building the tensor product of the one-dimensional rule at the finest level in every direction, Smolyak builds tensor products of fine levels in some directions and coarse levels in others, and combines many such products in a single rule. If the one-dimensional rule is nested, these tensor products have many coincident nodes – reducing the total cost. The resulting set of nodes is known as a sparse grid [19].

A concise formula for the Smolyak rule [26] is

[TABLE]

where $\{N_{k}\}_{k=1}^{N}\subset\mathbb{N}$ is an increasing sequence and $\mathcal{Q}_{N_{k}}$ is an $N_{k}$-node quadrature rule. $\{N_{k}\}$ is typically an exponentially growing sequence, because then the Smolyak cubature rule has a relatively high degree, which can be seen in the following lemma [18, 19].

Lemma 1**.**

Let $\{N_{k}\}$ grow exponentially in $k$. Then $S_{K}$ has at least degree $2(K-d)+1$.

In this paper, all Smolyak cubature rules are generated using quadrature rule sets with the following exponentially growing numbers of nodes:

[TABLE]

This sequence is chosen such that the sequence of $N_{k}$ nodes of the Clenshaw–Curtis quadrature rule is nested.

Smolyak rules do not have positive weights in general, but the condition number $\kappa$ of the cubature rule is bounded if the original quadrature rule has positive weights [19]:

[TABLE]

In Figure 1 two-dimensional sparse grids resulting from Smolyak applied to Gauss–Legendre, Gauss–Jacobi (with $\alpha=\beta=4$), and Clenshaw–Curtis quadrature rules are plotted. Thanks to nesting, Clenshaw–Curtis rules result in less than half of the number of nodes of the other rules, and this benefit will improve in higher dimensions. If Gaussian rules are used, the weights of the original rule are certainly positive, and therefore (3.5) holds.

4 Carathéodory reduction of quadrature rules

In this section new quadrature rules are introduced based on the removal of nodes. Given an initial quadrature rule with positive weights, a set of nested quadrature rules with positive weights is determined by removing nodes while retaining symmetry. We call these new rules reduced quadrature rules and the procedure to remove nodes the reduction step. These rules are by construction nested and do have positive weights. The reduction step can be applied in such a way that the rules are also symmetric.

4.1 Reduction step

The principle of the reduction step is as follows. First, recall the linear system (3.1):

[TABLE]

which describes an $N$-node quadrature rule of degree $N-1$. The goal is to find a subset of $N-1$ nodes that form a quadrature rule of degree $N-2$. Such a rule can easily be determined by considering the following system:

[TABLE]

where $V_{-1}$ is the matrix $V$ after the removal of the last row. This notation is used hereafter in a more general way: $A_{-k}$ denotes matrix $A$ after the removal of the last $k$ rows.

Each column of $V_{-1}$ is related to a node of the quadrature rule, so removing a column from the matrix above and solving the resulting system yields a nested quadrature rule of degree $N-2$. The question remains which column can be removed such that the system that remains has a solution with positive elements. The answer follows from (a variant of) the well-known Carathéodory theorem. The constructive proof will be useful later.

Theorem 1** (Carathéodory’s theorem).**

Let $\mathbf{v}_{1},\mathbf{v}_{2},\dots,\mathbf{v}_{N},\mathbf{v}_{N+1}$ be $N+1$ vectors spanning an $N$-dimensional space. Let $\mathbf{v}=\sum_{k=1}^{N+1}\lambda_{k}\mathbf{v}_{k}$ with $\lambda_{k}\geq 0$. Then there exist $\beta_{k}\geq 0$ such that $\mathbf{v}=\sum_{k\in I}\beta_{k}\mathbf{v}_{k}$ and $I\subset\{1,\dots,N+1\}$ with $|I|\leq N$.

Proof.

Because $\mathbf{v}_{1},\dots,\mathbf{v}_{N+1}$ are $N+1$ vectors in an $N$-dimensional space, they must be linearly dependent. So there are $c_{k}$, not all equal to zero, such that

[TABLE]

So for any $\alpha\in\mathbb{R}$, it is true that

[TABLE]

Without loss of generality, we assume that at least one $c_{k}>0$. Then the following choice is well-defined:

[TABLE]

Choosing $\beta_{k}=\lambda_{k}-\alpha c_{k}$, it is true that $\beta_{k_{0}}=0$ so with $I=\{1,2,\dots,k_{0}-1,k_{0}+1,\dots,N\}$ the following holds:

[TABLE]

Carathéodory’s theorem can be interpreted as a column-removal step. First, let $\mathbf{v}$ be the columns of $V_{-1}$. The elements $c_{k}$ from the proof form a null vector of the matrix. Determining $\alpha$ and $k_{0}$ from the proof yields that $w_{k}-\alpha c_{k}\geq 0$ and $w_{k_{0}}-\alpha c_{k_{0}}=0$, such that the node $x_{k_{0}}$ can be removed from the quadrature rule. This yields a quadrature rule of $N-1$ nodes of degree $N-2$ with positive weights, which was the goal.

Repeatedly applying the reduction step to an existing quadrature rule yields a set of nested quadrature rules with positive weights. The reduction step is however not unique in general. The null vector $\mathbf{c}$ contains both positive and negative elements (guaranteed by the fact that the first row of the matrix contains only positive values), so $-\mathbf{c}$ is also a null vector with both positive and negative elements and each null vector can be used to eliminate a different node.

This non-uniqueness imposes a choice. We suggest a heuristic greedy strategy of eliminating at each reduction step that node with the lowest probability based on the underlying probability density function (of the two nodes that can be eliminated). Nodes with high probability are retained. We call this the prior criterion. We shall see that for symmetric distributions and rules, this criterion does not apply (nodes have equal probability), which will be discussed in the next section.

The Smolyak cubature rule has been determined for several sets of reduced Gauss quadrature rules (see Figure 2, here the standard normal distribution is used) using the prior criterion. If two nodes have equal probability, the node which is most far from the center is removed. Comparing this to the Smolyak cubature rules which were determined previously (see Figure 1) yields that the number of nodes is the same as for the Clenshaw–Curtis quadrature rule, but the weights are positive (hence (3.5) can be used) and the location of the nodes is dependent on the distribution.

In the plots of the quadrature rules (below the sparse grids in Figure 2) it is clearly visible that the quadrature rules are not symmetric but do have positive weights. We extend the reduction step such that the rules are symmetric.

4.2 Symmetry

The reduction step does not necessarily keep a symmetric quadrature rule symmetric, because nodes are generally removed one-by-one. This is undesirable and does not happen if the null vector used in the procedure has the same symmetry as the weights, because eliminating one weight then automatically eliminates the symmetric weight.

Such a symmetric null vector always exists, as shown in the following lemma. The key notion is that to keep a symmetric quadrature rule symmetric two nodes have to be removed, which can be implemented by removing two rows from the Vandermonde-matrix (i.e. constructing $V_{-2}$ instead of $V_{-1}$). The proof is constructive, providing an algorithm for the reduction.

Lemma 2**.**

There exists a symmetric null vector of $V_{-2}$.

Proof.

The principle of the proof is to (i) construct a matrix $V^{\prime}_{-2}$ that encodes the symmetry property, (ii) prove that this matrix is singular and (iii) state a procedure to derive the null vector of this matrix.

Let $\{\xi_{k}\}_{k=1}^{N}$ be the nodes of an $N$-node symmetric quadrature rule.

A case distinction is made. First, let $N$ be even. Without loss of generality, we assume that $\xi_{1}<\xi_{2}<\dots<\xi_{N}$ and that the quadrature rule is symmetric around [math]. Then the nodes can be written as follows:

[TABLE]

Consider the following matrix:

[TABLE]

The matrix $V^{\prime}_{-2}$ is constructed by combining columns of $V_{-2}$. Each column of $V_{-2}$ is used exactly once. $V^{\prime}_{-2}$ is not a square matrix, so it is not trivial to see that a non-trivial null vector exists. However, if a null vector $\mathbf{c}^{\prime}$ exists, it can easily be transformed into a symmetric null vector of $V_{-2}$ using $\mathbf{c}={(c^{\prime}_{1},c^{\prime}_{2},\dots,c^{\prime}_{\frac{N}{2}},c^{\prime}_{\frac{N}{2}},\dots,c^{\prime}_{2},c^{\prime}_{1})}^{\mathrm{T}}$.

There always exists such a null vector $\mathbf{c}^{\prime}$ because for $p$ odd, it is true that

[TABLE]

Therefore, determining a null vector of $V^{\prime}_{-2}$ is equivalent to determining a null vector of the following matrix:

[TABLE]

Here, the rows of $V_{-2}$ consisting of zeros are removed and it is used that $\xi_{1}^{2}=\xi_{N}^{2}$, $\xi_{2}^{2}=\xi_{N-1}^{2}$, etc. $A_{-2}$ is an $\left(\frac{N}{2}-1\right)\times\frac{N}{2}$-matrix, which is singular, hence always has a non-trivial null vector.

If $N$ is odd, the same principle can be applied with the nodes:

[TABLE]

Therefore, after constructing $V^{\prime}_{-2}$ a null vector needs to be determined of the following matrix:

[TABLE]

In this case, $A_{-2}$ is an $\left\lfloor\frac{N}{2}\right\rfloor\times(\left\lfloor\frac{N}{2}\right\rfloor+1)$-matrix, which is again singular.

Concluding, in both cases ($N$ odd or even) there exists a symmetric null vector of $V_{-2}$, that can be constructed using a null vector of matrix $A_{-2}$. ∎

The previous lemma leads to the following main theorem about reduced symmetric quadrature rules with positive weights. Note that it is important to remove two nodes each time in both cases ($N$ odd or even), because only then it is possible to start with a quadrature rule of odd length and iteratively remove nodes until a quadrature rule of only the middle node is obtained.

Theorem 2**.**

Let $\{\xi_{1},\dots,\xi_{N}\}$ form an $N$-node symmetric quadrature rule with positive weights of degree $N-1$. Then there exist $\xi_{i}$ and $\xi_{j}$ with $i\neq j$ such that $\{\xi_{1},\dots,\xi_{N}\}\setminus\{\xi_{i},\xi_{j}\}$ forms an $(N-2)$-node symmetric quadrature rule with positive weights of degree $N-3$.

Proof.

The proof follows directly from Lemma 2. Let $A_{-2}$ as in Lemma 2 be given and let $\mathbf{c}$ be a null vector of $A_{-2}$. There are two cases:

• •

If $N$ is even: both $\mathbf{c}$ and $-\mathbf{c}$ yield the removal of two nodes in the reduction step.

• •

If $N$ is odd: it is possible that either $\mathbf{c}$ or $-\mathbf{c}$ yields the removal of the middle node, which would result in the removal of just one node. However, either $\mathbf{c}$ or $-\mathbf{c}$ yields the removal of two nodes.

In both cases, pick $i=k_{0}$ and $j=N-k_{0}+1$, where $k_{0}$ is from the proof of Carathéodory’s theorem. ∎

From this theorem and the case distinction between quadrature rules of even and odd length, an algorithm can be formulated that generates the nested quadrature rule keeping weights positive and a symmetric quadrature rule symmetric (see Algorithm 1). Although the matrices $V_{-1}$ and $A_{-2}$ can become ill-conditioned for large $N$, we did not observe any numerical issues in determining null vectors of these matrices for $N$ up to $2^{10}$.

Using this algorithm, symmetric reduced quadrature rules can be generated using the prior criterion. The resulting Smolyak cubature rules are therefore also symmetric (see Figure 3 for examples). The symmetry is also clearly visible in the plots of the quadrature rules.

5 Reduced cubature rules

In the previous section, a procedure has been outlined to generate a set of symmetric and nested quadrature rules with positive weights. Exactly the same principles can be applied to cubature rules, i.e. in a multi-dimensional setting.

The set-up is the same as in the previous section. First, the reduction step is introduced ignoring symmetry. This extension is straightforward. Then symmetries in multi-dimensional spaces are studied and a similar theory as in the one-dimensional case is developed regarding the symmetry of nested cubature rules.

5.1 Multi-dimensional reduction step

Let $\{\boldsymbol{\xi}_{1},\dots,\boldsymbol{\xi}_{N}\}$ and $\{w_{1},\dots,w_{N}\}$ be $N$ cubature nodes and positive weights in a $d$-dimensional space forming a cubature rule of degree $K$, with $N=\dim\mathbb{P}(K,d)$. Let $G$ be the $N\times N$ generalized Vandermonde-matrix, as introduced previously (see (3.2)). The goal is to determine a subset of nodes that forms a cubature rule of degree $K-1$, with positive weights. Such a cubature rule has a generalized Vandermonde-matrix of size $\dim\mathbb{P}(K-1,d)\times\dim\mathbb{P}(K-1,d)$.

The one-dimensional reduction step can be easily generalized to determine this nested cubature rule as follows. First, let $G$ be the generalized Vandermonde matrix again. Then $G_{-C}$ with $C=\dim\mathbb{P}(K,d)-\dim\mathbb{P}(K-1,d)$ has a $C$-dimensional null space†††Recall that $A_{-k}$ is matrix $A$ without its last $k$ rows.. Applying Carathéodory’s theorem iteratively to this matrix allows for the removal of $C$ columns, which yields the pursued $(\dim\mathbb{P}(K-1,d))$-node cubature rule with positive weights.

Just as in the one-dimensional case, the choice of basis for the null space is not unique, and since the number of null vectors is larger in multiple dimensions, more freedom to select the node to be removed is available. Iteratively applying the prior criterion is again an option, and will again result in loss of symmetry.

In multi-dimensional spaces, a symmetric $K$-degree cubature rule with positive weights of $\dim\mathbb{P}(K,d)$ nodes is not trivial to derive. However, a good initial cubature rule can be determined using the introduced reduction step. Starting with a $K$-degree tensor product rule, nodes can be removed from this rule until $\dim\mathbb{P}(K,d)$ nodes are left and using this cubature rule, a set of nested cubature rules can be generated.

5.2 Symmetries

As in the one-dimensional case, the reduction step does not keep a symmetric cubature rule symmetric. In a multi-dimensional space, many different types of symmetries can be considered. We consider two types of reflectional symmetry:

Symmetry along an axis, i.e. the plane of symmetry has the property $x^{(k)}=0$, where $x^{(k)}$ is a coordinate. If a cubature rule is symmetric in this way in all dimensions, the planes of symmetry divide the space into $2^{d}$ orthants (multi-dimensional quadrants). We call this a type-1 symmetry (see Figure 4 for a sketch). 2. 2.

Symmetry along a plane having $x^{(k)}=x^{(j)}$, where $x^{(k)}$ and $x^{(j)}$ are two coordinates. If a cubature rule is symmetric in this way in all dimensions, the planes of symmetry divide the space into $2^{d}$ orthants after a rotation over $\frac{1}{4}\pi$ of the complete basis. We call this a type-2 symmetry (see Figure 4 for a sketch).

In a tensor product cubature rule, the first symmetry occurs if the rule is generated using a symmetric quadrature rule. The second symmetry occurs if one quadrature rule is used multiple times in several dimensions.

To preserve symmetry after the removal of nodes, the null vector used to remove the nodes must have the same symmetry. As in the one-dimensional case, it is not guaranteed that such a null vector exists. Under certain conditions such a null vector does exist. The theory is more cumbersome than in the one-dimensional case, but has the same general structure: to determine a symmetric null vector of $G_{-C}$, a matrix $G^{\prime}$ is constructed and a proof is given that a null vector of $G^{\prime}$ can be transformed into a null vector of $G_{-C}$. Dependencies in the row space of $G^{\prime}$ finish the proof.

5.2.1 Type-1 symmetry

First we demonstrate the existence of a suitable null vector. The proof has the same structure as the proof of Lemma 2.

Lemma 3**.**

Let $\{\boldsymbol{\xi}_{1},\dots,\boldsymbol{\xi}_{N}\}$ be a type-1 symmetric cubature rule of degree $K$ with positive weights $\{w_{1},\dots,w_{N}\}$. Let $Q$ be an orthant and let $N_{\bar{Q}}$ be the number of cubature nodes in $\bar{Q}$. Then there exists a symmetric null vector of $G_{-C}$ if

[TABLE]

Proof.

See A. ∎

Using this lemma, a theorem can be stated about nested type-1 symmetric cubature rules with positive weights.

Theorem 3**.**

Let $\{\boldsymbol{\xi}_{1},\dots,\boldsymbol{\xi}_{N}\}$ be a type-1 symmetric cubature rule of degree $K$ with positive weights. Assume there are no cubature nodes shared between orthants (i.e., on the plane of symmetry). Then there exist $I=2^{d}\binom{\lfloor\frac{K-1}{2}\rfloor+d}{d}$ indices $i_{1},i_{2},\dots,i_{I}$ such that $\{\boldsymbol{\xi}_{i_{1}},\boldsymbol{\xi}_{i_{2}},\dots,\boldsymbol{\xi}_{i_{I}}\}$ forms a type-1 symmetric cubature rule of degree $K-1$ with positive weights.

Proof.

For a visual proof in two dimensions, see Figure 4. Let $N_{Q}$ be the number of nodes in an orthant $Q$. Because no nodes are shared between orthants, it is true that

[TABLE]

where $N_{\bar{Q}}$ is the number of nodes in $\bar{Q}$. Therefore the total number of nodes of the cubature rule equals $N_{Q}2^{d}$. From Lemma 3 it is known that the number of nodes that can be removed from one orthant equals

[TABLE]

Hence, the number of nodes remaining in the orthant equals

[TABLE]

So, the total number of nodes after all removal steps is $\binom{\lfloor\frac{K-1}{2}\rfloor+d}{d}2^{d}$. ∎

If there are cubature rule nodes shared between orthants (which is almost always the case), then $I$ is an upper bound of the number of nodes after a removal procedure.

5.2.2 Type-2 symmetry

Again, we state a lemma about the existence of a null vector. And again, the proof has the same structure as the proof of Lemma 2. However, the number of independent rows cannot be deduced explicitly anymore, such that the following lemma is necessary.

Lemma 4**.**

Let $\mathbf{s}=(s_{1},\dots,s_{d})\in\mathbb{N}_{+}^{d}$ be a sequence. If

• •

$\|\mathbf{s}\|_{1}\leq B$, where $B>0$ and $B\in\mathbb{N}$,

• •

$\mathbf{s}$* is weakly increasing, i.e. $s_{1}\leq s_{2}\leq\dots$*

*then there exist $1+\sum_{l=1}^{B}p_{d}(l)$ such sequences, where $p$ is the restricted partition function‡‡‡There are several definitions of the restricted partition number. Here, it is the number of compositions of the number $l$ with at most $d$ summands.. *

Lemma 5**.**

Let $\{\boldsymbol{\xi}_{1},\dots,\boldsymbol{\xi}_{N}\}$ be a type-2 symmetric cubature rule of degree $K$ with positive weights $\{w_{1},\dots,w_{N}\}$. Let $Q$ be an orthant after a rotation over $\frac{1}{4}\pi$ of all axes. Let $N_{\bar{Q}}$ be the number of cubature nodes in $\bar{Q}$. Then there exists a symmetric null vector of $G_{-C}$ if

[TABLE]

where $p_{d}(l)$ is the restricted partition function.

Proof.

See A. ∎

A similar theorem can be developed about the nested cubature rule in this case.

Theorem 4**.**

Let $\{\boldsymbol{\xi}_{1},\dots,\boldsymbol{\xi}_{N}\}$ be a type-2 symmetric cubature rule of degree $K$ with positive weights. Assume there are no cubature nodes shared between orthants after a rotation over $\frac{1}{4}\pi$ (i.e., on the plane of symmetry). Then there exist $I=2^{d}\left(1+\sum_{l=1}^{K-1}p_{d}(l)\right)$ indices $i_{1},i_{2},\dots,i_{I}$ such that $\{\boldsymbol{\xi}_{i_{1}},\boldsymbol{\xi}_{i_{2}},\dots,\boldsymbol{\xi}_{i_{I}}\}$ forms a type-2 symmetric cubature rule of degree $K-1$ with positive weights.

Proof.

Combine the proof of Theorem 3 with Lemma 5. For a visual proof in two dimensions, see Figure 4. ∎

If there are nodes shared between orthants, the theorem provides an upper bound of the number of nodes. The two lemmas can be combined into the following corollary. The resulting theorem has the same structure as the two theorems above and is therefore omitted.

Corollary 1**.**

Let $\{\boldsymbol{\xi}_{1},\dots,\boldsymbol{\xi}_{N}\}$ be a type-1 and type-2 symmetric cubature rule of degree $K$ with positive weights $\{w_{1},\dots,w_{N}\}$. Let $Q_{1}$ be an orthant and let $Q_{2}$ be an orthant after a rotation over $\frac{1}{4}\pi$ of all axes. Let $N_{\bar{Q_{1}}\cap\bar{Q_{2}}}\eqqcolon N_{\bar{Q}}$ be the number of nodes in both $\bar{Q_{1}}$ and $\bar{Q_{2}}$. Then there exists a symmetric null vector of $G_{-C}$ if

[TABLE]

where $p_{d}(l)$ is the restricted partition function.

Proof.

Combine Lemma 3 and 5. ∎

5.3 Reduced cubature rule

With the construction above, three different reduced cubature rules can be considered:

The reduced cubature rule, which is a set of nested cubature rules with positive weights (but no symmetry). This set can be generated by repeatedly applying the reduction step and is the multi-dimensional extension of the reduced quadrature rule. See Figure 5 for an example. 2. 2.

The symmetric reduced cubature rule, which is a set of nested cubature rules with positive weights, incorporating symmetry. This set can be generated by applying reduction steps with the theories deduced above. This is the multi-dimensional extension of the symmetric reduced quadrature rule. See Figure 5 for an example. 3. 3.

The negative symmetric reduced cubature rule, which is a set of nested cubature rules incorporating symmetry, but having possibly multiple negative weights. Although the weights are not absolutely bounded, this cubature rule yields a very small number of nodes. The set can be generated using the theories above trying to remove as many nodes as possible in each step. In the one-dimensional case, this rule was not studied because the number of nodes was not relevant. See Figure 5 for an example. In the figure it is clearly visible that nodes on the boundaries of the orthants of the theorems above are maintained.

The cubature rules are constructed using the multi-dimensional Vandermonde-matrix and incorporating Lemma 3 and 5. Pseudo-code for this is provided in B. The algorithm stated there considers a fixed fraction of a tensor grid, whose number of nodes increases rapidly. Therefore determining a null vector of the corresponding generalized Vandermonde-matrix becomes rapidly computationally expensive. The condition number of the matrix depends on the symmetries of the original cubature rule, i.e., the more symmetries there are, the better the condition number is. Creating a general efficient implementation is ongoing research.

If the number of dimensions is not too large ($d\lesssim 5$), then the cubature rule with positive weights has approximately the same number of nodes as the Smolyak sparse grid (see Figure 5 for an example). For higher dimensions, all cubature rules suffer from the curse of dimensionality. The cubature rule with positive weights has the largest growth in number of nodes. The cubature rule with (some) negative weights has the smallest growth compared with both the Smolyak cubature rule and the cubature rule with (only) positive weights.

Both the cubature rules with positive weights and negative weights are based on the removal of groups of nodes. If all groups would be of equal size, both approaches yield an equal number of nodes. However, nodes on the plane of symmetry belong to smaller groups. The cubature rule with positive weights removes nodes such that the weights remain positive and does not take this into account. The growth of the nodes with respect to the dimension is therefore large. In Figure 5 it is clearly visible that nodes on the plane of symmetry are being removed. As opposed to this, the cubature rule with negative weights does not remove these nodes (see Figure 5). Choosing a quadrature rule of odd length results therefore in a smaller cubature rule as there are more groups of nodes with equal weights.

Just as for the Smolyak cubature rule, the number of nodes cannot be deduced analytically but must be tabulated (see Table 1). For higher dimensional cases, the results for the positive reduced cubature rule are omitted due to computational constraints. In the table the differences in growth are clearly visible. In the 5-dimensional case it can be observed that choosing a quadrature rule of odd length yields less nodes. Selecting a larger initial rule can therefore result in a smaller reduced rule. For the same reason there are less results for the 15-, 20-, and 25-dimensional cases: the cases where the original quadrature rule has even length are computationally unfeasible.

5.4 Condition number

For the Smolyak cubature rule the growth of the condition number with respect to the number of nodes is bounded (recall (3.5)). The reduced quadrature rule has positive weights, so the condition number equals $1$ in this case. For the negative reduced cubature rule no such bounds exist, as far as the authors know. The condition number can be deduced numerically to assess its growth (see Figure 6). The maximum degree (which is 15 here) is chosen such that (numerically) the sum of the weights equals 1 with a maximum error of $10^{-12}$. We are primarily interested in $\kappa$, not in the numerical accuracy of the procedure.

The growth of the condition number of the Smolyak rules is equal, which is evident. The Smolyak rule generated with a reduced quadrature rule of the prior criterion has larger condition number than the Smolyak rule generated using Clenshaw–Curtis rules. If this is unwanted, we suggest a weight criterion where the reduced rule is selected with the smallest mutual difference, i.e., with the smallest $\max_{k}w_{k}-\min_{k}w_{k}$. The condition number of the Smolyak rule generated with this reduced quadrature rule is significantly smaller and close to the condition number of a Clenshaw–Curtis. However, this criterion does not use the underlying distribution, so we do not study it further.

The condition number of the reduced negative cubature rule is smaller than that of the Smolyak rule. However, for larger degrees severe numerical errors occur in the algorithm to generate these rules, which is not the case for the Smolyak rule.

6 Numerical results

In this section the proposed cubature rules are applied to multiple problems and compared with tensor product and Smolyak cubature rule.

This section is built as follows: in the first sub-section the cubature rules will be used to integrate the Genz test functions. These functions are designed for testing cubature rules. The second and third sub-section contain applications of the cubature rules to two UQ cases. In the second sub-section the standard lid-driven cavity flow problem with uncertain boundary conditions and material properties will be studied, using a Lattice Boltzmann method to compute the flow. In the final sub-section the main advantage of allowing negative weights is shown, i.e. high accuracy for a moderately high-dimensional problem. An aircraft aerodynamics test case is considered with seven uncertain parameters, using the Euler equations of gas dynamics and a finite-volume discretization of these to compute the corresponding aircraft aerodynamics. Because the conventional methods require a large number of simulations, only the results of the reduced cubature rule with negative weights are discussed in this case.

6.1 Genz test functions

6.1.1 Uniform distribution

To test the quality of cubature rules, several functions have been developed by Genz [8]. Each function has a different specific property or attribute, of which the effect can be enlarged by a parameter $\mathbf{a}$. A shape parameter $\mathbf{u}$ can be used to transform the function without changing the property (see Table 2 for all functions and their relevant attributes). For all functions the exact value of the integral can be determined [21].

Reducing a cubature rule only maintains the polynomial accuracy, which requires sufficient smoothness of the integrand. The first four Genz functions are in $C^{\infty}({[0,1]}^{d})$, while the fifth is in $C^{0}({[0,1]}^{d})$, and the sixth is only piecewise continuous. Hence a priori we expect spectral convergence for the first four functions and poor convergence for the fifth and sixth, independent of the particular rule.

To obtain meaningful, instructive results the coefficients $\mathbf{a}$ and $\mathbf{u}$ are chosen randomly, with each component from similar uniform distributions, subject to the constraints $\|\mathbf{a}\|_{2}=2.5$ and $\|\mathbf{u}\|_{2}=1$. Moreover, each component of both $\mathbf{a}$ and $\mathbf{u}$ is positive. The integration error is determined with respect to the exact solution, and averaged over 100 runs.

Convergence plots for all methods and all Genz functions are depicted in Figure 7. The multi-dimensional reduced rules are initiated using tensor products of Gaussian quadrature rules and the reduction procedure is only applied once to keep numerical artifacts small. The Smolyak procedure is applied twice using Clenshaw–Curtis quadrature rules or reduced quadrature rules using a fine Clenshaw–Curtis rule as initial rule. The results can be divided into three classes that exhibit different behaviors: (i) $f_{1}$, $f_{2}$, $f_{4}$, (ii) $f_{3}$, and (iii) $f_{5}$, $f_{6}$.

Class (i) is formed by smooth results, which show almost spectral convergence for all methods. Both symmetric reduced rules consistently outperform the tensor product, and negative symmetric reduced rules also consistently outperform both Smolyak rules. No significant difference exists between the two Smolyak rules.

Class (ii) is an exception, most likely caused by the concentration of mass at one corner of the integration domain. Both Smolyak and reduced (negative weights) rules remove nodes at corners, and thereby poorly approximate the most important region of the integrand. The reduced quadrature rule keeps some nodes at the corner up to small levels, so therefore the Smolyak rule with reduced rules performs slightly better. The tensor and reduced (positive weights) rules do have nodes there, which makes the error much smaller.

For class (iii) the integrands lack sufficient smoothness for polynomial approximations to be stable. As expected, spectral convergence is not evident, but some limited linear convergence is visible. In both cases Smolyak rules acquit themselves well compared to all other methods.

In summary the proposed reduced rules are empirically converging at the level of Smolyak or better, given sufficient smoothness in the integrand.

6.1.2 Non-uniform distribution

The reduced cubature rule can be determined for any distribution whose moments can be evaluated. Therefore we assess the convergence of the rules using a $\beta(10,10)$-distribution, a highly non-uniform distribution. The Smolyak cubature rule with Clenshaw–Curtis quadrature rules is not considered anymore. We only examine the results of $f_{1}$, $f_{2}$, and $f_{4}$, i.e., the cases where convergence was observed.

The results are created in a similar way as in the uniform case, i.e., the coefficients $\mathbf{a}$ and $\mathbf{u}$ are chosen randomly subject to the constraints $\|\mathbf{a}\|_{2}=2.5$ and $\|\mathbf{u}\|_{2}=1$. The integration error is determined with respect to a reference value, calculated using a $30^{5}$ tensor grid created with Gaussian quadrature rules. Convergence plots are depicted in Figure 8.

In comparison with the previous results we see better convergence of $f_{1}$ for all cubature rules. This is due to the $\beta(10,10)$-distribution, which damps the oscillations of the function under consideration. Moreover it is clearly visible that an approximate value of the integral is used here.

$f_{2}$ and $f_{4}$ show similar results: in both cases the reduced rule with negative weights shows the best results. The differences are larger in this case, which is due to the $\beta(10,10)$-distribution that introduces many small weights in the initial quadrature and cubature rule. These weights are prone to removal in both the negative and positive reduction algorithm.

6.1.3 Dimension dependence

All integrals so far have been determined using 5-dimensional rules. To quantify the performance of rules depending on the dimension, the integration error is studied for varying dimension.

The results are again created in a similar way as in the previous case, i.e. using the random coefficients and averaging the result. The uniform distribution is reconsidered, such that an exact value of the integral is known. We again limit ourselves to $f_{1}$, $f_{2}$, and $f_{4}$. All cubature rules are generated such that they are of degree 9. The number of nodes is not taken into account here (but can be found in Table 1). The integration errors up to 10 dimensions are plotted in Figure 9.

The oscillatory function $f_{1}$ and the Gaussian function do not become more difficult to integrate in higher dimensions, as the Taylor expansions are comparable to the one-dimensional case. However, the mass of the product peak of $f_{2}$ becomes smaller as the dimension increases, which makes the integral easier to evaluate numerically, as integrating the peak accurately becomes of less importance. This is also reflected in the results. The tensor product rule shows excellent results for both functions, which is due to the exact integration of more polynomials in comparison to the other rules.

The positive reduced cubature rule shows the smallest growth compared to the other rules (excluding the tensor rule). This is due to the positive weights, that yield a condition number equal to 1.

The negative reduced cubature rule and the two Smolyak rules show similar growth. Although the negative reduced rule does not reduce or increase the integration error, it does need much less nodes to obtain this error compared with both Smolyak rules (see Table 1).

In summary, the tensor grid yields the largest cubature rule and shows the smallest error. The reduced positive rule needs less nodes, but yields a larger error. The two Smolyak rules and the reduced negative rule have approximately equivalent error, but the reduced negative rule yields a much smaller grid.

6.2 Lid-driven cavity flow test case, using a Lattice Boltzmann method

6.2.1 Problem description

The standard lid-driven cavity flow (e.g. [11]) is considered with two uncertain flow parameters, and four UQ methods are compared, namely MC, SC with a Smolyak sparse grid, and SC with the two new reduced cubature rules. A sketch of the geometry and the imposed boundary conditions is given in Figure 10. The boundary condition imposed at both singular corners is $\mathbf{u}=0$, where $\mathbf{u}$ is the fluid velocity vector.

The deterministic problem is solved using a Lattice Boltzmann method. The implementation used for the current case is a straightforward D2Q9 BGK-model using Zou–He boundary conditions [33]. Reference data for several values of the Reynolds number can be found in Ghia et al. [11]. The results from the Lattice Boltzmann implementation compare well with the data provided (see Figure 11).

Two uncertain parameters are specified (see Table 3). Both parameters have a $\beta(a,b)$-distribution, with probability density function

[TABLE]

The ranges are chosen such that the Reynolds number based on the lid velocity is between 10 and 400 (see Figure 12 for the solutions of the two extreme cases).

6.2.2 Results

To evaluate the accuracy of the methods we consider the $u$-component of the fluid velocity everywhere in the domain. A reference mean solution is obtained using a fine tensor product rule of $65\times 65$ Gaussian nodes in the parameter space, resulting in a reference mean solution $\bar{u}^{*}$, depicted in Figure 13 (the tensor grid is plotted in Figure 14).

UQ is applied with four different methods:

MC using random samples; 2. 2.

SC with a Smolyak sparse grid created with positive symmetric reduced Gauss–Jacobi quadrature rules, see Figure 14; 3. 3.

SC with the symmetric reduced rule initiated with a tensor grid with positive weights, see Figure 14; 4. 4.

SC with the negative symmetric reduced rule, see Figure 14.

The number of nodes is chosen in such a way that the degree of the resulting rule equals 13. For the Smolyak cubature rule, this can be achieved by choosing $K=8$, because then $2(K-d)+1=13$ (see Lemma 1). The initial quadrature rule is chosen such that it is the finest quadrature rule used by the Smolyak procedure. For the reduced cubature rules, the initial cubature rule is a $13\times 13$ tensor grid of Gaussian quadrature rules. All grids (except the MC nodes) are shown in Figure 14.

For a measure of accuracy we take the $L^{2}$-norm of the difference between the predicted mean velocity field $\bar{u}(N)$ and the reference mean field $\bar{u}^{*}$, where $N$ is the number of nodes, i.e.

[TABLE]

The convergence is shown in Figure 15. The $\mathcal{O}(1/\sqrt{N})$ convergence of MC is already becoming apparent. The polynomial-based methods seem to show spectral convergence, suggesting that the response is smooth with respect to the parameters, as might be expected from physical considerations. Of the polynomial methods, the Smolyak rule and symmetric reduced rule perform approximately the same. The negative symmetric reduced rule performs poorly. The authors attribute this to the low dimension of the problem; we have seen the benefits of allowing negative weights primarily in five and more dimensions. In low-dimensional cases the difference in the number of nodes is too small to overcome the large absolute differences in weights (see the color bar in Figure 14).

6.3 Aircraft aerodynamics test case, using the Euler equations and a finite-volume method

6.3.1 Problem description

To show the performance the reduced cubature rule with negative weights (i.e. high accuracy with a relatively small number of nodes) we consider an aircraft aerodynamics test case with seven uncertain parameters. Only the reduced cubature rule with negative weights is used as the other cubature rules require too many nodes.

The geometry of the airplane is based on a sample airplane geometry of the program sumo [6], the so-called “twin-engine utility aircraft” (see Figure 16).

As flow model we consider the Euler equations of gas dynamics. The Euler-flow problem is solved using the second-order accurate finite-volume code $\text{SU}^{2}$ [20]. The tools sumo and TetGen are used for mesh generation [23]. Besides modeling the surfaces, sumo creates surface meshes, which are used as input for TetGen which generates the volume meshes with a spherical far field boundary. See Figure 17 for an example solution.

Of the seven uncertain parameters, three are geometrical and assumed to be normally distributed. The mean is the base geometry value and the standard deviation is defined to be 5%. Four uncertain operational parameters are considered in addition, all modeled as $\beta(4,4)$ variables. See Table 4 for details.

The geometrical uncertain parameters are specifically chosen such that the 4-digit NACA airfoil series can be used to parameterize them. The three parameters of the NACA series are essentially these parameters. The base geometry, which defines the mean of the distributions of these parameters, is that of the NACA2412 airfoil.

6.3.2 Results

A reduced tensor cubature rule with negative weights of degree 9 is generated, which yields a cubature rule of 1,293 nodes. A Smolyak sparse grid of the same degree consists of 2,465 nodes. If a tensor grid is used, then 78,125 simulations are necessary (if Gaussian rules are used). The reduced cubature rule with positive weights consists of 8,713 nodes in this case. This example shows that if time is an issue, allowing negative weights can indeed reduce the number of nodes significantly.

The lift, drag, and side-force coefficients are scalars obtained by integration over the wetted surface. Their moments are listed in Table 5 together with estimations based on least-squares regression. The degree of the polynomial fitted using least-squares is 5, which is the maximum number possible to keep the system determined, i.e.

[TABLE]

Although the cubature rule has negative weights the variance is non-negative and the values are close to the least-squares estimates. Moreover, the order of magnitude of the lower-order moments seems physically reasonable (although no reference data is available for this case).

We calculate moments of the pressure coefficient over the entire surface, and plot them in Figure 18. Here, the four $\beta$-distributed uncertain parameters are taken into account, because geometrical uncertainties cannot be plotted. The mean looks very much like a typical pressure distribution. The variance (which is non-negative at all nodes on the geometry) highlights the location of the shocks, for which a small change in location leads to a large change in the pressure distribution. This result is consistent with UQ analyses of airfoil flows with shocks [29]. The higher-order moments (which have shown convergence due to the removal of the three geometrical uncertainties) are also large near the shock, indicating that it is unlikely that the resulting distribution is Gaussian.

7 Conclusion

Non-intrusive uncertainty quantification was studied using stochastic collocation methods. Three important properties are relevant for quadrature rules and cubature rules if they are used in stochastic collocation: nested, positive weights, and symmetry if the original distribution is symmetric. Existing quadrature rules and cubature rules do not have these three properties. The Clenshaw–Curtis quadrature rule is symmetric and nested, but does not have positive weights in general. The Gauss quadrature rule is symmetric and has positive weights, but is not nested.

The introduced quadrature rule performs well using a Smolyak grid. Any positive quadrature rule can be used to generate a set of nested quadrature rules, which can be used as input for a Smolyak procedure. Gaussian rules are quadrature rules which are always positive.

The proposed cubature rule comes in three variants: one which ignores symmetry, one which ignores positive weights, and one that has all properties. For low-dimensional problems, the cubature rule which satisfies all three properties has approximately the same nodes as a Smolyak sparse grid, but has positive weights. Convergence is also approximately equal, which can be seen in the results from the Genz test functions and the lid-driven cavity flow. If the uncertainty quantification problem is high-dimensional and computational efficiency is important, the positivity constraint can be relaxed to remove more nodes. The symmetric reduced cubature rule with possibly some negative weights yields less nodes than Smolyak grids and the reduced cubature rule with positive weights. In our example, higher-order moments remained positive, although the cubature rule has negative weights. This can become an issue though in applications where the response surface is more complex.

Although the cubature rules are nested, it is also important to note that the initial cubature rule still influences the result. If the initial cubature rule is too small, no general strategy exists to add nodes to the cubature rule and create a larger one, having the three properties. This is an option for further research.

Appendix A Proofs of Lemmas 3 and 5

A.1 Proof of Lemma 3

Lemma**.**

Let $\{\boldsymbol{\xi}_{1},\dots,\boldsymbol{\xi}_{N}\}$ be a type-1 symmetric cubature rule of degree $K$ with positive weights $\{w_{1},\dots,w_{N}\}$. Let $Q$ be an orthant and let $N_{\bar{Q}}$ be the number of cubature nodes in $\bar{Q}$. Then there exists a symmetric null vector of $G_{-C}$ if

[TABLE]

Proof.

Without loss of generality, assume that the cubature rule is symmetric around 0. First we construct a suitable matrix $G^{\prime}$.

Let $Q$ be an orthant. Without loss of generality, let $\{\boldsymbol{\xi}_{1},\dots,\boldsymbol{\xi}_{N_{\bar{Q}}}\}\in\bar{Q}$. Let $G^{\prime}$ be the generalized Vandermonde-matrix of the nodes

[TABLE]

omitting the monomials with an odd power. $G^{\prime}$ is a $\binom{\lfloor\frac{K}{2}\rfloor+d}{d}\times N_{\bar{Q}}$-matrix.

Assume that

[TABLE]

holds. Let $\mathbf{c}^{\prime}$ be a null vector of $G^{\prime}$. Construct vector $\mathbf{c}$ as follows for $k=1,\dots,N$:

[TABLE]

where $j$ is such that

[TABLE]

and $1\leq j\leq N_{\bar{Q}}$. Due to the symmetry, such a $j$ always exists. Now $\mathbf{c}$ is a null vector of $G$. To see this, let a row index $i$ of $G$ be given, with row $\mathbf{r}_{i}$. A case distinction is made.

Case 1: $m_{i}$ only contains even powers. Let $m_{i}(\boldsymbol{\xi})=\boldsymbol{\xi}^{\alpha}$. Then

[TABLE]

because this monomial was included in matrix $G^{\prime}$.

Case 2: $m_{i}$ contains an odd power. Without loss of generality, assume that the first power is odd, i.e., $m_{i}(\boldsymbol{\xi})=\boldsymbol{\xi}^{\alpha}$, with $\alpha^{(1)}$ odd. Then let $I$ be the index set of cubature nodes with first element equal to 0, $J$ the index set of cubature nodes with first element larger than 0, and $K$ all other indices. Due to type-1 symmetry, the size of $J$ and $K$ is equal. Split the nodes along the plane of symmetry. Then

[TABLE]

A.2 Proof of Lemma 5

Lemma 6**.**

Let $\{\boldsymbol{\xi}_{1},\dots,\boldsymbol{\xi}_{N}\}$ be a type-2 symmetric cubature rule of degree $K$ with positive weights $\{w_{1},\dots,w_{N}\}$. Let $Q$ be an orthant after a rotation over $\frac{1}{4}\pi$ of all axes. Let $N_{\bar{Q}}$ be the number of cubature nodes in $\bar{Q}$. Then there exists a symmetric null vector of $G_{-C}$ if

[TABLE]

where $p_{d}(l)$ is the restricted partition function.

Proof.

Without loss of generality, assume that the cubature rule is symmetric around 0. The proof has the same structure as the proof of Lemma 2, combined with Lemma 6. Again, we first construct a suitable matrix $G^{\prime}$.

Let $Q$ be an orthant after $\frac{1}{4}\pi$ rotation of all basis vectors. Without loss of generality, let $\{\boldsymbol{\xi}_{1},\dots,\boldsymbol{\xi}_{N_{\bar{Q}}}\}\in\bar{Q}$. Let $G^{\prime}$ be the generalized Vandermonde-matrix of the nodes

[TABLE]

where $\#\sigma\boldsymbol{\xi}_{k}$ is the number of permutations of the elements of cubature node $\boldsymbol{\xi}_{k}$ and where $\sigma_{k}$ is the $k^{\text{th}}$ permutation operator, i.e., it is a sum over all permutations of $\boldsymbol{\xi}_{k}$. Omit all monomials $\boldsymbol{\xi}^{\alpha}$ with a power which is not sorted. Due to Lemma 4, $G^{\prime}$ is a $\left(1+\sum_{l=1}^{K}p_{d}(l)\right)\times N_{\bar{Q}}$-matrix. Assume that

[TABLE]

holds. Let $\mathbf{c}^{\prime}$ be a null vector of $G^{\prime}$. Construct vector $\mathbf{c}$ as follows for $k=1,\dots,N$:

[TABLE]

where $j$ is such that

[TABLE]

for a suitable permutation $\sigma$ with $1\leq j\leq N_{\bar{Q}}$. This is well-defined because of the type-2 symmetry. Now $\mathbf{c}$ is a null vector of $G$. To see this, let a row index $i$ of $G$ be given, with row $\mathbf{r}_{i}$. Let $m_{i}(\boldsymbol{\xi})=\boldsymbol{\xi}^{\alpha}$ be the respective monomial. Then:

[TABLE]

For each $\boldsymbol{\xi}_{k}$, there exists a permutation operator $\sigma_{k}$, such that $\sigma_{k}(\boldsymbol{\xi}_{k})=\boldsymbol{\xi}_{j}$ for $1\leq j\leq N_{\bar{Q}}$. Hence:

[TABLE]

Appendix B Algorithms for reduced cubature rule generation

The algorithms provided in this appendix are not yet efficient for high-dimensional cubature rules of high degree, and are merely given for sake of completeness. Examples of possible optimizations are the following: numerical issues can arise when determining the null vectors and time and memory issues can arise if the complete $G_{-C}$ is constructed. The first issue can be partially overcome by scaling all elements of the nodes onto the same interval and the second issue can be circumvented by constructing the matrix $G_{-C}$ column-wise and checking for existence of a null vector after each addition of a column.