Numerical Integration in Multiple Dimensions with Designed Quadrature
NNUMERICAL INTEGRATION IN MULTIPLE DIMENSIONS WITHDESIGNED QUADRATURE ∗ VAHID KESHAVARZZADEH † , ROBERT M. KIRBY † , ‡ , AND
AKIL NARAYAN † , § Abstract.
We present a systematic computational framework for generating positive quadra-ture rules in multiple dimensions on general geometries. A direct moment-matching formulationthat enforces exact integration on polynomial subspaces yields nonlinear conditions and geometricconstraints on nodes and weights. We use penalty methods to address the geometric constraints,and subsequently solve a quadratic minimization problem via the Gauss-Newton method. Our anal-ysis provides guidance on requisite sizes of quadrature rules for a given polynomial subspace, andfurnishes useful user-end stability bounds on error in the quadrature rule in the case when the poly-nomial moment conditions are violated by a small amount due to, e.g., finite precision limitationsor stagnation of the optimization procedure. We present several numerical examples investigatingoptimal low-degree quadrature rules, Lebesgue constants, and 100-dimensional quadrature. Our cap-stone examples compare our quadrature approach to popular alternatives, such as sparse grids andquasi-Monte Carlo methods, for problems in linear elasticity and topology optimization.
Key words.
Numerical Integration, Multi Dimensions, Polynomial Approximation, QuadratureOptimization
AMS subject classifications.
1. Introduction.
Numerical quadrature, the process of computing approxima-tions to integrals, is widely used in many fields of science and engineering. A conve-nient and popular choice is a quadrature rule that uses point evaluations of a function f : (cid:90) Γ f ( x ) ω ( x )d x ≈ n (cid:88) j =1 f ( x j ) w j , where Γ is some set in d -dimensional Euclidean space R d , ω is a positive weightfunction, and x j and w j are the nodes and weights, respectively, of the quadraturerule that must be determined. The main desirable properties of quadrature rulesare accuracy for a broad class of functions, a small number n of nodes/weights, andpositivity of the weights. (Positive weights are desired so that the absolute conditionnumber of the quadrature rule is controlled.)In one dimension, Gaussian quadrature rules [29, 44] satisfy many of these desir-able properties, but computing an efficient quadrature rule (or “cubature” rule) forhigher dimensions is a considerably more challenging problem. When Γ and ω areof tensor-product form, one straightforward construction results from tensorization ∗ Accepted for publication in SIAM Journal on Scientific Computing - Methods and Algorithmsfor Scientific Computing section.
Funding:
This research was sponsored by ARL under Cooperative Agreement NumberW911NF-12-2-0023. The views and conclusions contained in this document are those of the au-thors and should not be interpreted as representing the official policies, either expressed or implied,of ARL or the U.S. Government. The U.S. Government is authorized to reproduce and distributereprints for Government purposes notwithstanding any copyright notation herein. The first andthird authors are partially supported by AFOSR FA9550-15-1-0467. The third author is partiallysupported by DARPA EQUiPS N660011524053. † Scientific Computing and Imaging Institute, University of Utah, Salt Lake City, UT ([email protected], [email protected], [email protected]). ‡ School of Computing, University of Utah, Salt Lake City, UT § Department of Mathematics, University of Utah, Salt Lake City, UT1
This manuscript is for review purposes only. a r X i v : . [ c s . NA ] A p r V. KESHAVARZZADEH, R. M. KIRBY, AND A. NARAYAN of univariate quadrature rules. However, the computational complexity required toevaluate f at the nodes of a tensorized quadrature rule quickly succumbs to the curseof dimensionality.Substantial progress has been made in constructing attractive multivariate quadra-ture rules. Sparse grids rely on a sophisticated manipulation of univariate quadra-ture rules [7, 16]. Quasi-Monte Carlo methods generate sequences that have low-discrepancy properties [33, 34, 37]. Mathematical characterizations of quadraturerules with specified exactness on polynomial spaces yield efficient nodes and weights[5, 9, 42, 56].The main contribution of this paper is a systematic computational approach fordesigning multivariate quadrature rules with exactness on general finite-dimensionalpolynomial spaces. Using polynomial exactness as a desideratum for constructingquadrature rules is not the only approach one could use (e.g., quasi-Monte Carlomethods do not adopt this approach). However, when the integrand f can be accu-rately approximated by a polynomial expansion with a small number of significantterms, then approximating the integral with a quadrature rule that is designed tointegrate the significant terms can be very efficient [10, 11]. In particular, finite-dimensional polynomial spaces can well-approximate solutions to some parametricoperator equations [12], and empirical tests with many engineering problems showthat polynomial approximations are very efficient [1, 2, 8].Our computational approach revolves around optimization; many algorithms forcomputing nodal sets via optimization have already been proposed [29, 30, 36, 45, 46,48, 53]. Our method, which we call designed quadrature , has the following advantages: • we can successfully compute nodal sets in up to 100 dimensions; • positivity of the weights is ensured; • quadrature rules over non-standard geometries can be computed; and • a prescribed polynomial accuracy can be sought over general polynomialspaces, not restricted to, e.g., total degree spaces.Our approach is simple: we formulate moment-matching conditions and geometricconstraints that prescribe nonlinear conditions on the nodes and weights. This directformulation allows significant flexibility with respect to geometry, weight function ω , and polynomial accuracy. Indeed, our procedures can compute quadrature ruleswith hyperbolic cross polynomial spaces, see Section 4.5, and can constrain nodallocations to awkward geometries, see Section 4.4. Our computational approach isto use constrained optimization algorithms to compute a quadrature rule from themoment-matching conditions. Our mathematical analysis provides a stability boundon error of the quadrature rule if the moment-matching conditions are violated (e.g.,due to numerical finite precision). We apply our designed quadrature rules to severalrealistic problems in computational science, including problems in linear elasticityand topology optimization. Comparisons against competing methods, such as sparsegrids and low-discrepancy sequences, illustrate that designed quadrature often attainssuperior accuracy with many fewer nodes.Our procedure is not without shortcomings: Being a direct moment-matchingproblem, our framework relies on large-scale optimization in high dimensions. For aspecified polynomial subspace on which we require integration accuracy, we cannot apriori determine the number of nodes that our procedure will produce (although wereview some theory that provides upper and lower bounds for n ). We likewise cannotensure that our algorithm produces an optimal quadrature rule size, but our numericalresults suggest favorable comparison with alterative techniques, see Section 4.2. Someof the optimization tools we use have tunable parameters; we have made automated This manuscript is for review purposes only.
ESIGNED QUADRATURE
2. Multivariate Quadrature.2.1. Notation.
Let ω be a given non-negative weight function (e.g., a probabilitydensity function) whose support is Γ ⊂ R d , where d ≥ x ∈ R d has components x = (cid:0) x (1) , x (2) , . . . , x ( d ) (cid:1) . The space L ω (Γ) is the setof functions f defined by L ω (Γ) = (cid:8) f : Γ → R (cid:12)(cid:12) (cid:107) f (cid:107) < ∞ (cid:9) , (cid:107) f (cid:107) = ( f, f ) , ( f, g ) = (cid:90) Γ f ( x ) g ( x ) ω ( x )d x . We use standard multi-index notation: α ∈ N d denotes a multi-index, and Λ acollection of multi-indices. We have α = ( α , . . . , α d ) , x α = d (cid:89) j =1 (cid:16) x ( j ) (cid:17) α j , | α | = d (cid:88) j =1 α j . We impose a partial ordering on multi-indices via component-wise comparisons: with α , β ∈ N d , then α ≤ β if and only if all component-wise inequalities are true. Amulti-index set Λ is called downward closed if α ∈ Λ = ⇒ β ∈ Λ ∀ β ≤ α . We assume throughout this paper that the weight function has finite polynomialmoments of all orders: (cid:90) Γ ( x α ) ω ( x ) < ∞ , α ∈ N d . This assumption ensures existence of polynomial moments. Our ultimate goal is toconstruct a set of n points { x q } nq =1 ⊂ Γ and positive weights w q > I ( f ) = (cid:90) Γ f ( x ) ω ( x )d x ≈ n (cid:88) q =1 w q f ( x q ) , (1a)for functions f within a “large” class of functions. We attempt to achieve this byenforcing equality above for f in a subspace Π of polynomials: (cid:90) Γ f ( x ) ω ( x )d x = n (cid:88) q =1 w q f ( x q ) , f ∈ Π . (1b)The quadrature strategy is accurate if f can be well-approximated by a polynomialfrom Π. There are numerous technical conditions on Π and f that yield quantitativestatements about polynomial approximation accuracy, e.g., [3]. In this article, weassume that Π is given and fixed through some a priori study ensuring that there This manuscript is for review purposes only.
V. KESHAVARZZADEH, R. M. KIRBY, AND A. NARAYAN exists a polynomial in Π that accurately approximates f to within some user-specifiedtolerance. Typically we will define Π through some finite multi-index set Λ:Π = span (cid:8) x α (cid:12)(cid:12) α ∈ Λ (cid:9) . In many applications, the function f typically exhibits smoothness (e.g., inte-grable high-order derivatives), which in turn implies that polynomial approximationsconverge at a high order with respect to the degree of approximation. Under the as-sumption that f is smooth, we therefore expect that the integral of a polynomial thatapproximates f to be a good approximation if the approximating polynomial spaceΠ contains high-degree polynomials. Our main goal in this paper is then familiarwhen viewed through the lens of classical analysis: make Π as large as possible whilekeeping n as small as possible.Two particularly popular choices for polynomial spaces Π can be defined by theindex setsΛ T r = (cid:8) α ∈ N d (cid:12)(cid:12) | α | ≤ r (cid:9) , Λ H r = α ∈ N d (cid:12)(cid:12) d (cid:89) j =1 ( α j + 1) ≤ r + 1 , for some non-negative integer r . Both of these multi-index sets are downward closed.The total order and hyperbolic cross polynomial subspaces are defined by, respectively,Π T r = span (cid:8) x α (cid:12)(cid:12) α ∈ Λ T r (cid:9) , Π H r = span (cid:8) x α (cid:12)(cid:12) α ∈ Λ H r (cid:9) . (2) The algorithm we present in this paper applies to general polynomial spaces, butour numerical examples will focus on the spaces above since they are common inlarge-scale computing problems. When Γ ⊂ R , the optimal quadra-ture rule is provided by the ω -Gauss quadrature rule. In one dimension, we use theshorthand Π k = Π T k . The first step in defining this rule is to prescribe an orthonor-mal basis for Π k . A Gram-Schmidt argument implies that such a basis of orthonormalpolynomials exists with elements p m ( · ), where deg p m = m . All univariate orthonor-mal polynomial families satisfy the three-term recurrence relation, xp m ( x ) = (cid:112) b m p m − ( x ) + a m p m ( x ) + (cid:112) b m +1 p m +1 ( x ) , (3)for m ≥
0, with p − ≡ p ≡ / √ b to seed the recurrence. The recurrencecoefficients are given by a m = ( xp m , p m ) , b m = ( p m , p m )( p m − , p m − ) , for m ≥
0, with b = ( p , p ). Classical orthogonal polynomial families, such as theLegendre and Hermite polynomials, fit this mold with explicit formula for the a n and b n coefficients [44]. Gaussian quadrature rules are n -point rules that exactly integratepolynomials in Π n − [39, 14]. Theorem
Let x , . . . , x n be the roots of the n th or-thogonal polynomial p n ( x ) and let w , . . . , w n be the solution of the system of equations (4) n (cid:88) q =1 p j ( x q ) w q = (cid:40) √ b , if j = 00 , for j = 1 , . . . , n − . This manuscript is for review purposes only.
ESIGNED QUADRATURE Then x q ∈ Γ and w q > for q = 1 , , . . . , n and (5) (cid:90) Γ ω ( x ) p ( x ) dx = n (cid:88) q =1 p ( x q ) w q holds for all polynomials p ∈ Π n − . Historically significant algorithmic strategies for computing Gauss quadrature rulesare given in [15, 18]. The elegant linear algebraic formulations described in thesereferences compute the quadrature rule with knowledge of only of a finite number ofrecurrence coefficients a n , b n . If Γ and ω ( x ) are both tensorial, then thegeneralization of univariate orthogonal polynomials to multivariate ones is straight-forward. The tensorial structure impliesΓ = × dj =1 Γ j , ω ( x ) = d (cid:89) j =1 ω j (cid:16) x ( j ) (cid:17) , for univariate domains Γ j ⊂ R and univariate weights ω j ( · ). If p ( j ) n ( · ) is the univariateorthonormal polynomial family associated with ω j over Γ j , then π α ( x ) = d (cid:89) j =1 p ( j ) α j (cid:16) x ( j ) (cid:17) , α ∈ N d , (6)defines a family of multivariate polynomials orthonormal under ω , i.e., ( π α , π β ) = δ α , β , where δ is the Kronecker delta. The polynomial spaces in (2) can be written asΠ T r = span (cid:8) π α (cid:12)(cid:12) α ∈ Λ T r (cid:9) Π H r = span (cid:8) π α (cid:12)(cid:12) α ∈ Λ H r (cid:9) The following result is the cornerstone of our algorithm:
Proposition
Let Λ be a multi-index set with ∈ Λ . Suppose that x , . . . , x n and w , . . . , w n are the solution of the system of equations (7) n (cid:88) q =1 π α ( x q ) w q = (cid:40) /π , if α = , if α ∈ Λ \{ } then (8) (cid:90) Γ ω ( x ) π ( x ) d x = n (cid:88) q =1 π ( x q ) w q holds for all polynomials π ∈ Π Λ . The proof is straightforward by noting that (cid:82) Γ π α ( x ) ω ( x )d x = 0 when α (cid:54) = dueto orthogonality, and thus (7) is a moment-matching condition. Unlike Theorem 2.1,this multivariate result does not guarantee the positivity of weights nor does it ensurethat the nodes lie in Γ. We enforce these conditions in our computational frameworkin Section 3. Finally, we note that Proposition 2.2 is true even when Γ and ω are nottensorial. We concentrate on the tensorial situation in this paper because a tensorialassumption is standard for large dimension d . This manuscript is for review purposes only.
V. KESHAVARZZADEH, R. M. KIRBY, AND A. NARAYAN
One of the main uses of quadrature rules is in the construction of polynomialapproximation via discrete quadrature. If f is a given continuous function and Θ is agiven multi-index set, then f ( x ) ≈ f Θ ( x ) = (cid:88) α ∈ Θ (cid:98) f α π α ( x ) , (cid:98) f α = n (cid:88) q =1 π α ( x q ) f ( x q ) w q , (9)where (cid:98) f α are meant to approximate the Fourier ( L ω -projection) coefficients of f .Ideally, if f ∈ Π Θ then f Θ = f , i.e., this construction reproduces polynomials in Π Θ .As one expects, this only happens when the quadrature rule is sufficiently accurate,as defined by the size of Λ in (7). Proposition
Let Λ be a downward-closed multi-index set, and suppose that x q and w q for q = 1 , . . . , n define a quadrature rule satisfying (7) . Let Θ be any indexset satisfying Θ + Θ = (cid:8) α + β (cid:12)(cid:12) α , β ∈ Θ (cid:9) ⊆ Λ . (10) If f ∈ Π Θ , then f Θ defined in (9) satisfies f Θ = f .Proof. Suppose f ∈ Π Θ , so that f ( x ) = (cid:88) α ∈ Θ f α π α ( x ) , f α = ( f, π α ) , where the formula for the coefficients f α is due to orthogonality. We will show thatthe computed quadrature coefficients (cid:98) f α defined in (9) satisfy (cid:98) f α = f α . Fix β ∈ Θ.Then, f ( x ) π β ( x ) = (cid:88) α ∈ Θ f α π α ( x ) π β ( x ) . There are coefficients c α , γ such that π α = (cid:88) γ ≤ α c α , γ x γ . Therefore, π α ( x ) π β ( x ) = (cid:88) γ ≤ α c α , γ x γ (cid:88) γ ≤ β c β , γ x γ = (cid:88) γ ≤ α + β d α , β , γ x γ , for some coefficients d α , β , γ . The index α + β ∈ Λ owing to the assumption (10), andsince Λ is downward closed, then we have that π α ( x ) π β ( x ) ∈ Π Λ . Therefore, the n -point quadrature rule integrates π α ( x ) π β ( x ), and thus (cid:98) f β = n (cid:88) q =1 f ( x q ) π β ( x q ) = (cid:88) α ∈ Θ f α n (cid:88) q =1 π α ( x ) π β ( x ) = (cid:88) α ∈ Θ f α ( π α , π β ) = f β , Since (cid:98) f β = f β , then f Θ = f . This manuscript is for review purposes only.
ESIGNED QUADRATURE n -point Gauss quadrature rule,we can reproduce polynomials up to degree n −
1: Take Λ = { , . . . , n − } , andchoose Θ = { , . . . , n − } . The polynomial f Θ constructed by the procedure (9)matches the function f if f ∈ Π Θ since Θ + Θ ⊂ Λ. The above result codifies thiscondition in the multivariate case. Note that Θ ⊂ Λ is not a strict enough conditionsince the approximate Fourier coefficients defined in (9) will not necessarily be accu-rate. We also note that the integrand is a product of polynomials, therefore requiringexactness on polynomial products is the correct condition, hence the Θ + Θ ⊂ Γrequirement.Given a multi-index set Λ, there is a smallest possible quadrature size n such that(7) holds. This smallest n is given by the size of the largest Θ satisfying (10). Theorem
Let Λ be a downward-closed index set. The size n of anyquadrature rule satisfying (7) has lower bound n ≥ L (Λ) := max (cid:8) | Θ | (cid:12)(cid:12) Θ + Θ ⊆ Λ (cid:9) . The number L (Λ) defined above is called the maximal half-set size in [24], and a cor-responding L (Λ)-point quadrature rule is a minimal rule. In that reference, concreteexamples of (i) non-existence, and of (ii) existence but non-uniqueness of minimal mul-tivariate quadrature rules achieving the lower bound above are shown. If Λ = Λ n − in the univariate case, Gaussian quadrature rules are non-unique. Our numerical al-gorithm essentially seeks to find minimal rules, but we can rarely find such quadraturerules. However, our generated quadrature rule sizes are only modestly larger than theoptimal L (Λ). Gaussian quadrature rules defined by Theorem2.1 can be computed via linear algebra, but multivariate quadrature rules definedby (7) have no known analogous computational simplification. In order to solve thisnonlinear system of equations we utilize Newton’s method. We therefore expect that(7) is not exactly satisfied by the computed solution, or it is satisfied to within sometolerance.Fixing a downward-closed index set Λ with size M = | Λ | , consider the matrix X ∈ R d × n whose n columns are the samples x j , and let w ∈ R n be a vector containingthe n weights. Let V ( X ) ∈ R n × M denote the Vandermonde-like matrix with entries( V ) k,j = π α ( k ) ( x j ) , j = 1 , . . . , n, k = 1 , . . . , M, (11)where we have introduced an ordering α (1) , . . . α ( m ) on the elements of Λ. We assume α (1) = , but the remaining ordering of elements is irrelevant. The system (7) canthen be written as V ( X ) w = e /π , where e = (1 , , , . . . , T ∈ R M is a cardinal unit vector. Instead of achieving theequality above, our computational solver computes an approximate solution ( X , w )to the above system, satisfying(12) (cid:107) V ( X ) w − e /π (cid:107) = (cid:15) ≥ . Our next result quantifies the effect of the residual (cid:15) on the accuracy of the designedquadrature rule. To prove this result, we require the additional assumption that thequadrature weights are positive, which is enforced in our computations.
This manuscript is for review purposes only.
V. KESHAVARZZADEH, R. M. KIRBY, AND A. NARAYAN
Proposition
Let ω ( x ) be a probability density function on Γ , and let Λ beany multi-index set containing (i.e., Π Λ contains constant functions). Assume that ( X , w ) satisfies (12) with some (cid:15) ≥ , and assume the weights are all positive. Thenfor any f ∈ L ω (Γ) , (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:90) f ( x ) ω ( x )d x − n (cid:88) q =1 w q f ( x q ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ≤ (cid:15) (cid:107) f (cid:107) + max j =1 ,...n | f ( x j ) − p ( x j ) | , (13) where p ∈ Π Λ is the L ω (Γ) -orthogonal projection of f onto Π Λ . This result does apply to all our computed designed quadrature rules since we enforcepositivity of the weights. It is not applicable to other polynomial-based rules whereweights can be negative, such as sparse grids.
Proof.
For an arbitrary p ∈ Π Λ , the following holds p ( x ) = (cid:88) α ∈I p α π α ( x ) , p α = ( p, π α ) , (14)and thus (cid:107) p (cid:107) = ( p, p ) = (cid:80) α ∈ Λ p α . We have: (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:90) Γ f ( x ) ω ( x )d x − n (cid:88) q =1 w q f ( x q ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ≤ (cid:12)(cid:12)(cid:12)(cid:12)(cid:90) Γ ( f ( x ) − p ( x )) ω ( x )d x (cid:12)(cid:12)(cid:12)(cid:12)(cid:124) (cid:123)(cid:122) (cid:125) ( a ) + (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) n (cid:88) q =1 ( p ( x q ) − f ( x q )) w q (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:124) (cid:123)(cid:122) (cid:125) ( b ) (15) + (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:90) Γ p ( x ) ω ( x )d x − n (cid:88) q =1 w q p ( x q ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:124) (cid:123)(cid:122) (cid:125) (c) We now choose p as the L ω (Γ)-orthogonal projection of f into Π Λ : p = argmin q ∈ Π Λ (cid:107) f − q (cid:107) = ⇒ (cid:90) Γ [ f ( x ) − p ( x )] φ ( x ) ω ( x )d x = 0 ∀ φ ∈ Π Λ . (16)Since ∈ Λ, the above holds in particular for φ ( x ) ≡ (cid:12)(cid:12)(cid:12)(cid:12)(cid:90) Γ ( f ( x ) − p ( x )) ω ( x )d x (cid:12)(cid:12)(cid:12)(cid:12) = 0Term (b) can be bounded as(b) ≤ n (cid:88) q =1 | w q | | p ( x q ) − f ( x q ) | ≤ max q =1 ,...,n | p ( x q ) − f ( x q ) | , where the last inequality uses the fact that (cid:80) Nq =1 | w q | = (cid:80) Nq =1 w q = (cid:82) Γ ω ( x )d x = 1since the weights are positive and ω is a probability density. Finally, term (c) can bebounded as follows: Since p ∈ Π Λ then by (14), n (cid:88) q =1 w q p ( x q ) = n (cid:88) q =1 (cid:88) α ∈ Λ w q p α π α ( x q ) = (cid:88) α ∈ Λ p α (cid:32) n (cid:88) q =1 w q π α ( x q ) (cid:33) This manuscript is for review purposes only.
ESIGNED QUADRATURE V ( X ) w fromthe relation (12); note also that (cid:98) π α cf. Equation (9) equals an entry in the vector b .Therefore, combining the above equation and using the Cauchy-Schwarz inequality:(c) = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:90) Γ p ( x ) ω ( x )d x − n (cid:88) q =1 w q p ( x q ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:88) α ∈ Λ p α (cid:32)(cid:90) Γ π α ( x ) ω ( x )d x − n (cid:88) q =1 w q π α ( x q ) (cid:33)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:88) α ∈ Λ p α (cid:32) δ α , /π − n (cid:88) q =1 w q π α ( x q ) (cid:33)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ≤ (cid:115) (cid:88) α ∈ Λ p α (cid:107) V ( X ) w − e /π (cid:107) ≤ (cid:15) (cid:107) p (cid:107) ≤ (cid:15) (cid:107) f (cid:107) , where the final inequality is Bessel’s inequality, which holds since we have chosen p as in (16). Combining our estimates for terms (a), (b), and (c) in (15) completes theproof.Relative to the pointwise error committed by best L ω (Γ) approximations, theestimate provided by Proposition 2.5 bounds the quadrature error in terms of thequantity (cid:15) , which is explicitly computable given a quadrature rule. A (Smolyak) sparse grid is astructured point configuration in multiple dimensions, formed from unions of ten-sorized univariate rules. Quadrature weights often accompany points in a sparse grid.We briefly describe sparse grids for polynomial integration in this section; they willbe used for comparison in our numerical results section.Consider a tensorial Γ as in Section 2.3, and for simplicity assume that the uni-variate domains Γ j = Γ are the same, and that the univariate weights ω j = ω arethe same. Let X i denote a univariate quadrature rule (nodes and weights) of “level” i ≥
1, and define X = ∅ . The number of points n i in the quadrature rule X i is in-creasing with i , but can be freely chosen. For multi-index i ∈ N d , a d -variate tensorialrule and its corresponding weights are(17) A d, i = X i ⊗ . . . ⊗ X i d , w ( q ) = d (cid:89) r =1 w ( q r ) i r The univariate difference operator between sequential levels is written as∆ i = X i − X i − , i ≥ , (18)and for any k ∈ N , this approximation difference can be used to construct a d -variate,level- k -accurate sparse grid operator [7, 38], A d,k = k − (cid:88) r =0 (cid:88) i ∈ N d | i | = d + r ∆ i ⊗ . . . ⊗ ∆ i d = k − (cid:88) r = k − d ( − k − − r (cid:18) d − k − − r (cid:19) (cid:88) i ∈ N d | i | = d + r X i ⊗ . . . ⊗ X i d , (19)where the latter equality is shown in [52]. If the univariate quadrature rule X i exactlyintegrates univariate polynomials of order 2 i − A d,k isexact for d -variate polynomials of total order 2 k − X i to obtain optimal efficiency, but since the differences ∆ i appear in the Smolyak construction, then instead utilizing nested univariate rules cangenerate sparse grids with many fewer nodes than non-nested constructions. One can This manuscript is for review purposes only. V. KESHAVARZZADEH, R. M. KIRBY, AND A. NARAYAN use, for example, nested Clenshaw-Curtis rules [55], the nested Gauss-Patterson orGauss-Kronrod rules [16, 26, 35], or Leja sequences [31].Sparse grids have been used with great success in many modern applications, andthus are a good candidate for comparison against our approach of designed quadrature.However, sparse grids that integrate polynomials in a certain multi-index set use farmore points than the minimum number prescribed by Theorem 2.4 (see Figure 3for an empirical comparison), and frequently produce quadrature rules with negativeweights. Our results in Section 4 show that designed quadrature uses many fewerpoints than sparse grids for a given accuracy level, and guarantees positive quadratureweights.
3. Computational Framework.
Our procedure aims to compute nodes X = { x , . . . , x n } ∈ Γ n and positive weights w ∈ (0 , ∞ ) n that enforce equality in (7). Adirect formulation of (7) is(20) R ( d ) = V ( X ) w − e /π = , x j ∈ Γ , j = 1 , . . . , nw j > , j = 1 , . . . n where d = ( X , w ) are the decision variables. Instead of directly solving this con-strained root finding problem, we introduce a closely related constrained optimizationproblem(21) min X , w || R || subject to x j ∈ Γ , j = 1 , . . . , nw j > , j = 1 , . . . , n Clearly a solution to (20) also solves (21), but the reverse is not necessarily true.We compute solutions to (21), and when these solutions exhibit large nonzero valuesof (cid:107) R (cid:107) , we increase the quadrature rule size n and repeat. Using this strategy, weempirically find that for a specified (cid:15) we can satisfy (cid:107) R (cid:107) ≤ (cid:15) in all situations we havetried. Thus, our approach solves a relaxed version of (20) via repeated applicationsof (21). Our computational approach to solve (21) requires four major ingredients,each of which are described in the subsequent sections:Section 3.1 – Penalization: objective augmentation, transforming constrained rootfinding into unconstrained minimization problemSection 3.2 – Iteration: unconstrained minimization via the Gauss-Newton algo-rithmSection 3.3 – Regularization: numerical regularization to address ill-conditionedGauss-Newton update stepsSection 3.4 – Initialization: specification of an initial guessWe highlight above that regularization is required for our optimization. The objective R in (21) is highly ill-conditioned as a function of the decision variables. Withoutregularization, the update steps specified by the Gauss-Newton algorithm generallydo not result in convergence. However, with the regularization, we have found thatour optimization results in steps with decreasing residual. These observations can becorroborated by the numerical results in Section 4, and in particular Table 2 that listsCPU time and iterations required for computing 4-dimensional rules.Since our algorithm only minimizes the norm of R , the quadrature rule we com-pute is not guaranteed to integrate any polynomials exactly, only up to some tolerance This manuscript is for review purposes only.
ESIGNED QUADRATURE (cid:15) ≥ (cid:107) R (cid:107) . This is the utility of Proposition 2.5: if our optimization algo-rithm terminates with a particular value of (cid:15) , we have a quantitative understandingof how (cid:15) affects the quality of the quadrature rule relative to best L -approximatingpolynomials.Since we produce a quadrature rule that is only (cid:15) -exact, there may be manyquadrature rules that achieve this tolerance. In particular, our algorithm is not guar-anteed to produce optimal quadrature rules, but in comparison with some other tab-ulated rules from [40, 43, 53, 54], we find that our nodal counts are no greater than inthose references. There is one lone exception for integrating degree-8 polynomials inthree dimensions, where we find a rule with one point greater than reported in [53].Details are in Section 4.2 and in Table 1.Finally, our algorithm is subject to the same limitations as many other minimiza-tion algorithms: it may only find a local minimum of the objective, and not a globalminimum. Penalty methods are techniques to solve constrained opti-mization problems such as (21). Penalty methods augment the objective with a highcost for constraint violated, and subsequently solve an unconstrained optimizationproblem on the augmented objective.We use a popular penalty function, the non-negative and smooth quadratic func-tion. For example in d = 1 dimensions on Γ = [ − ,
1] with an n -point quadraturerule, the constraints and corresponding penalties P j , j = 1 , . . . , ( d + 1) n = 2 n as afunction of the 2 n decision variables d = ( X , w ) can be expressed as − ≤ x j ≤ ⇒ P j ( d ) = (max[0 , x j − , − − x j ]) ,w j ≥ ⇒ P n + j ( d ) = (max[0 , − w j ]) , for j = 1 , . . . , n . The total penalty associated with the constraints is then P ( d ) = ( d +1) n (cid:88) j =1 P j ( d ) . A penalty function approach to solve the constrained problem (21) uses a sequenceof unconstrained problems indexed by k ∈ N having objective functions g ( c k , d ) := (cid:13)(cid:13)(cid:13) (cid:101) R k (cid:13)(cid:13)(cid:13) = (cid:107) R (cid:107) + c k P ( d ) , (22)where we have defined the vector (cid:101) R k = R c k P c k P ... c k P ( d +1) n . The positive constants c k are monotonically increasing with k , i.e., c k +1 > c k . Eachunconstrained optimization yields an updated solution point d k , and as c k → ∞ the solution point of the unconstrained problem will converge to the solution of con-strained problem. The following lemma, adopted from [28], is used to show conver-gence of the penalty method. This manuscript is for review purposes only. V. KESHAVARZZADEH, R. M. KIRBY, AND A. NARAYAN
Lemma
Let d k be the minimizer for g ( c k , · ) and c k +1 > c k . Then: g ( c k , d k ) ≤ g ( c k +1 , d k +1 ) , P ( d k ) ≥ P ( d k +1 ) , || R ( d k ) || ≤ || R ( d k +1 ) || . Furthermore, let d ∗ be a solution to problem (21) . Then for each k , || R ( d k ) || ≤ g ( c k , d k ) ≤ || R ( d ∗ ) || . The above lemma denotes that the sequence of g ( c k , d k ) is nondecreasing andbounded above by the optimal objective value of the constrained optimization prob-lem. The following theorem establishes the global convergence of the penalty method.More precisely it verifies that any limit point of the sequence is a solution to (21). Theorem
Let { d k } , k ∈ N be a sequence of minimizers of (22) .Then any limit point of the sequence is a solution to problem (21) i.e. lim k ∈ N P ( d k ) =0 and lim k ∈ N || R ( d k ) || ≤ || R ( d ∗ ) || . The above theorem shows both that a limit point denoted by ¯ d is a feasible solutionsince P ( ¯ d ) = 0, and that it is optimal since || R ( ¯ d ) || ≤ || R ( d ∗ ) || .We can now formulate an unconstrained minimization problem with sequence ofincreasing c k on the objectives g in (22) for the decision variables d = ( X , w ),min d g ( c k , d )(23)which replace the constrained root-finding problem (20).It remains for us to specify how the constants c k are chosen: if d is the currentiterate for the decision variables, we use the formula c k = max (cid:26) A, || R ( d ) || (cid:27) , where A is a tunable parameter that is meant to be large. We use A = 10 in oursimulations. Also note that we never have c k = ∞ so that our iterations cannot exactlyconstrain the computed solution to lie in the feasible set. To address this in practicewe reformulate constraints to have non-zero penalty within a small radius inside thefeasible set. For example, instead of enforcing w j >
0, we enforce w j > − .Note that one may also consider barrier/interior point methods to enforce con-straints; however, in our algorithm we find that penalty methods are more suitable intransforming the constrained root finding problem to an unconstrained minimizationproblem. Having transformed the constrained prob-lem (21) into a sequence of unconstrained problems (23), we can now use standardunconstrained optimization tools.Two popular approaches for unconstrained optimization are gradient descent andNewton’s method. Both approaches in the context of our minimization require theJacobian of the objective function with respect to the decision variables. We define (cid:101) J k = ∂ (cid:101) R k ∂ d = J c k ∂P /∂ d c k ∂P /∂ d ... c k ∂P ( d +1) n /∂ d , J ( d ) := ∂ R ∂ d ∈ R M × ( d +1) n , (24) This manuscript is for review purposes only.
ESIGNED QUADRATURE ∂P j ∂ d ∈ R × ( d +1) n is the Jacobian of P j with respect to the decision variables.With use of our quadratic penalty function, these penalty Jacobians are Lipschitzcontinuous in the decision variables, and easily evaluated since they are quadraticfunctions. The matrix J has entries( J ) m, ( i − d + j = ∂π α ( m ) ( x i ) ∂x ( j ) i w i , ( J ) m,nd + i = π α ( m ) ( x j ) , (25)for m = 1 , . . . , M , i = 1 , . . . , n , and j = 1 , . . . , d . Above, we define π α ( m ) as in (11).Computing entries of the Jacobian matrix J is straightforward: Assuming the basis π α is of tensor-product form, (see Section 2.3) then we need only compute derivativesof univariate polynomials. A manipulation of the three-term recurrence relation (3)yields the recurrence (cid:112) b m +1 p (cid:48) m +1 ( x ) = ( x − a m ) p (cid:48) m ( x ) − √ b m p (cid:48) m − ( x ) + p m ( x ) . The partial derivatives in J may be evaluated using the relation above along with (6).We index iterations with k , which is the same k as that defining the sequenceof unconstrained problems (23). Thus, our choice of c k changes at each iteration.Gradient descent proceeds via iteration of the form d k +1 = d k − α ∂ (cid:107) (cid:101) R k (cid:107) ∂ d , ∂ (cid:107) (cid:101) R k (cid:107) ∂ d = (cid:101) J Tk (cid:101) R (cid:107) (cid:101) R k (cid:107) , with α a customizable step length that is frequently optimized via, e.g., a line-searchalgorithm. In contrast, a variant of Newton’s root finding method applied to rectan-gular systems is the Gauss-Newton method [39], having update iteration d k +1 = d k − ∆ d , ∆ d = (cid:16) (cid:101) J Tk (cid:101) J k (cid:17) − (cid:101) J Tk (cid:101) R k , (26)where both (cid:101) J k and (cid:101) R k are evaluated at d k . The iteration above reduces to thestandard Newton’s method when the system is square, i.e., M = n ( d + 1). Newton’smethod converges quadratically to a local solution for a sufficiently close initial guess d versus the gradient descent which has linear convergence [4]. We find that Gauss-Newton iterations are robust for our problem.Assuming an initial guess d is given, we can repeatedly apply the Gauss-Newtoniteration (26) until a stopping criterion is met. We terminate our iterations when theresidual norm falls below a user-defined threshold (cid:15) i.e. || (cid:101) R || < (cid:15) .A useful quantity to monitor during the iteration process is the magnitude of theNewton decrement, which often reflects quantitative proximity to the optimal point[6]. In its original form, the Newton decrement is the norm of the Newton step inthe quadratic norm defined by the Hessian. I.e., for optimizing f ( x ), the Newtondecrement norm is || ∆ d || ∇ f ( x ) = (∆ d T ∇ f ( x )∆ d ) / , where ∇ f is the Hessian of f . In our minimization procedure with non-squared systems we use(27) η = (cid:0) ∆ d T ( (cid:101) J Tk (cid:101) R k ) (cid:1) / . as a surrogate for a Hessian-based Newton decrement which decreases as d → d ∗ .Finally we note that, for a given quadrature rule size n , we cannot guarantee thata solution to (20) exists. In this case our Gauss-Newton iterations will exhibit residual This manuscript is for review purposes only. V. KESHAVARZZADEH, R. M. KIRBY, AND A. NARAYAN norms stagnating at some positive value while the Newton decrement is almost zero.When this occurs, we re-initialize the decision variables and enrich the current set ofdecision variables with additional nodes and weights and continue the optimizationprocedure. This procedure of gradually increasing the number of nodes and weightsis described more in Section 3.4.
The critical part of our minimization scheme is the eval-uation of Newton step (26). For our rectangular system, this is the least-squaressolution ∆ d to the linear system (cid:101) J ∆ d = (cid:101) R , where (cid:101) J = (cid:101) J k (cid:0) d k (cid:1) , and (cid:101) R = (cid:101) R k (cid:0) d k (cid:1) ; in this section we omit explicit notationaldependence on the iteration index k . The matrix (cid:101) J is frequently ill-conditioned,which hinders a direct solve of the above least-squares problem. To address this wecan consider a generic regularization of the above equality:(28) minimize ∆ d || (cid:101) J ∆ d − (cid:101) R || p subject to || ∆ d || q < τ, where p , q , and τ are free parameters. The trade off between the objective norm andsolution norm is characterized as a Pareto curve and shown to be convex in [50, 51]for generic norms 1 ≤ ( p, q ) ≤ ∞ . Exploiting this Pareto curve, the authors in [50, 51]devise an efficient algorithm and implementation [49] for computing the regularizedsolution when p = 2 , q = 1. These values correspond to the LASSO problem [47],which promotes solution sparsity and subset selection.Since sparsity is not our explicit goal, we opt for p = q = 2. This problem can besolved exactly [19], but at significant expense and the procedure lacks clear guidanceon choosing τ . We thus adopt an alternative approach. A penalized version of the p = q = 2 optimization (28) is Tikhonov regularization:(29) ∆ d λ = argmin (cid:110) || (cid:101) J ∆ d − (cid:101) R || + λ || ∆ d || (cid:111) , where λ is a regularization parameter that may be chosen by the user. This parameterhas significant impact on the quality of the solution with respect to the original least-squares problem. Assuming that we have a definitive value for λ , then the solution to(29) can be obtained via the singular value decomposition (SVD) of (cid:101) J . The SVD ofmatrix (cid:101) J N × M (for N < M) is given by(30) (cid:101) J = N (cid:88) i =1 u i σ i v Ti . where σ i are singular values (in decreasing order), and u i and v k are the correspondingleft- and right-singular vectors, respectively. The solution ∆ d λ is then obtained as(31) ∆ d λ = N (cid:88) i =1 ρ i u Ti (cid:101) R σ i v i . where ρ i are Tikhonov filter factors denoted by(32) ρ i = σ i σ i + λ (cid:39) (cid:40) σ i (cid:29) λσ i /λ σ i (cid:28) λ This manuscript is for review purposes only.
ESIGNED QUADRATURE λ . Therefore a suitable λ is bounded by the extremal singular values of (cid:101) J . Oneapproach to select λ is via analysis of the “ L -curve” of singular values [20, 21]. Thecorner of L -curve can be interpreted as the point with maximum curvature; evaluationor approximation of the curvature with respect to singular value index can be usedto find the index with maximum curvature, and the singular value corresponding tothis index prescribes λ .In practice, we evaluate the curvature of the singular value spectrum via finitedifferences on log( σ i ) (where the singular values are directly computed) and selectthe singular value that corresponds to the first spike in the spectrum. The regular-ization parameter can be updated after several, e.g., 30, Gauss-Newton iterations.However, for small size problems, i.e., small dimension d and | Λ | , a fixed appropriate λ throughout the Gauss-Newton scheme also yields solutions.Based on our numerical observations, adding a regularization parameter to all sin-gular values and computing the regularized Newton step as ∆ d λ = (cid:80) N i =1 [( u Ti (cid:101) R ) / ( σ i + λ )] v i enhances the convergence when d is close to the root i.e. || (cid:101) R || is small. The first step of the algorithm requires an initial guess d for nodes and weights; a particularly difficult aspect of this is the initial choice ofquadrature rule size n . Our algorithm tests several values of quadrature rule sizes n between an upper and lower bound; the determination of these bounds are describedbelow.With the multi-index set Λ given, Theorem 2.4 provides a lower bound on thevalue of n , and this lower bound L (Λ) is the optimal size for a quadrature rule.We are unaware of sufficient conditions under which optimal quadrature rules exist.However, optimal-size quadrature rules have been shown in special cases, e.g., [41],for total degree spaces Λ T k with k = 2 , ,
5. We have found that our algorithm is ableto recover these optimal-sized rules in the previously-mentioned cases.We formulate an upper bound on quadrature rule sizes based on a popular com-petitor: sparse grid constructions. The number of sparse grid points | A d,k | requiredto satisfy (7) with Λ = Λ T k can be estimated as | A d,k | ≈ (2 d ) k − ( k − [13] for sparse gridconstructions with non-nested univariate Gauss quadrature rules. Tabulation of theexact number of points for sparse grids constructed via univariate nested rules fromthe Hermite and Legendre systems is provided in [22].Our numerical results show that the number of designed quadrature nodes neededto satisfy (7) is n = κ | A d,k | where κ ∈ [0 . , .
9] using | A d,k | from [22]. We have foundthat an effective approach to choose the number of points is to perform a backtrackingline-search procedure, which initializes κ = 0 .
9, solves the optimization problem, andgradually decreases κ until the Gauss-Newton method does not converge to a desirabletolerance. Our strategy for eliminating nodes when κ is decreased is to discard thosewith the smallest weights.After the initial pass that generates n nodes and weights achieving (cid:107) (cid:101) R (cid:107) ≤ (cid:15) , weattempt to remove nodes with smallest weights as described previously. However,this may cause the optimization to stagnate without achieving the desired tolerance.When this happens, we enrich the nodal set by gradually adding more nodes untilwe can achieve the tolerance. This process is repeated until the elimination andenrichment procedures result in no change of the quadrature rule size; see Algorithm1, lines 9-18.Once an initial number of nodes n is determined ( κ = 0 . d - This manuscript is for review purposes only. V. KESHAVARZZADEH, R. M. KIRBY, AND A. NARAYAN variate Monte Carlo samples or Latin Hypercube samples are generated as the initialnodes. This is easily done for the domain Γ = [ − , d . Weights can be generateduniformly at random [0 ,
1] with (cid:80) i w i = | Λ | or set as a fixed value, e.g., w i = | Λ | /n .We normalize the weights by | Λ | in the numerical procedure to avoid very smallweights. To accommodate for this, we can set 1 /π = | Λ | in (20), and after we obtaina solution we can re-normalize the weights based on the true value of 1 /π .On the domain Γ = R d , we are usually concerned with the weight ω ( x ) =exp( −(cid:107) x (cid:107) ). Monte Carlo samples can be generated as realizations of a standardnormal random variable, and we transform Latin Hypercube samples on [0 , d to R d via inverse transform sampling corresponding to a standard normal random variable.(When Λ contains polynomials of very high degree there are more sophisticated sam-pling methods that can produce better initial guesses [32].) We initialize the weightsby setting w i = exp( −|| x i || /
2) and normalizing w i with respect to | Λ | as describedabove.Algorithm 1 summarizes Sections 3.1–3.4, including all the steps for our designedquadrature method. Algorithm 1
Designed Quadrature Initialize nodes and weights d with n = 0 . | A d,k | and specify the residual toler-ance, e.g., (cid:15) = 10 − . Set n = 0. while || (cid:101) R || > (cid:15) do Compute (cid:101) R and (cid:101) J using (22), (20), (24), and (25). Determine the regularization parameter λ from the SVD of (cid:101) J Compute the regularized Newton step ∆ d from (31) Update the decision variables d k +1 = d k − ∆ d . Compute the residual norm || (cid:101) R || and Newton decrement η from (27). if η < (cid:15) and || (cid:101) R || (cid:29) (cid:15) then Increase n , initialize new nodes and weights, and go to line 3. end if end while if n = n then Return else n ← n Decrease n by eliminating nodes with smallest weights, go to line 3. (Seediscussion about κ in Section 3.4.) end if4. Numerical Examples.4.1. Illustrative numerical example in d=2. In this example we consider d = 2 for a uniform weight on Γ = [0 , with an r = 2 total degree polynomial spacewith index set Λ T . This index set has six indices, corresponding to six constraintsin (7). Using n = 3 nodes there are ( d + 1) n = 9 decision variables. Note that exactformulas for the optimal quadrature rule is known in this case [43]. The augmentedJacobian (cid:101) J in (24) is a 15 × , and use uniform weights. The singular values of the Jacobian matrixare shown in Figure 1 for the initial and final decision variables corresponding to This manuscript is for review purposes only.
ESIGNED QUADRATURE λ . The results suggest thatany positive value in [0 . ,
5] can be used as a regularization parameter. We use aconstant λ throughout the iterations and fix the residual tolerance (cid:15) = 10 − . Theevolution of residual (cid:107) (cid:101) R (cid:107) , Newton decrement η , and penalty parameter c k is shownin Figure 1. Smaller λ values appear to yield faster convergence. i σ i InitialFinal i σ i InitialFinal i σ i InitialFinal k -20 M agn i t ude || e R || η c k -10 || e R || η c k -10 || e R || η c Fig. 1 . Singular values of Jacobian (cid:101) J cf. Equation (24) (top) and residual norm (cid:107) (cid:101) R (cid:107) ,Newton decrement η and penalty parameter c with respect to iterations k (bottom) for regu-larization parameter λ = 0 . (left), λ = 1 (middle) and λ = 5 (right). To visualize the optimal points for this quadrature, we randomize the initialnode positions and compute designed quadrature for 100 initializations. Plots ofthe ensemble of converged quadrature rules in two and three dimensions are shownin Figure 2. A set of 3-points (initial and final design) in each experiment forms atriangle i.e. vertices of each triangle are the quadrature points where each triangle isvisualized for distinguishment. The cumulative time for 100 designs took ∼ sec withMATLAB on a single core personal desktop, and each design takes ∼
15 iterationswith λ = 1. In this example we considerthe number of nodes required to achieve exact polynomial accuracy on total degreespaces Λ T r of various orders and dimensions. Our goal is to compare designed quadra-ture against sparse grids. The number of nodes required for exact integration on asparse grid is from [22]. Our tests fix dimension d = 3 and sweep values of the order r , and fix r = 5 and sweep values of the dimension d . We present the nodal countsin Table 1 and in Figure 3. Table 1 shows that designed quadrature consistentlyresults in fewer nodes than sparse grids for moderate values of r and d . We againemphasize that the weights for designed quadrature are all positive, unlike sparse gridquadrature.Figure 3 compares various node counts: The number of nodes in the product ruleis simply n d where n is the number of univariate Gauss quadrature nodes and the“lower bound” is the value L (Λ) determined from Theorem (2.4). Using Theorem 2.1 This manuscript is for review purposes only. V. KESHAVARZZADEH, R. M. KIRBY, AND A. NARAYAN x (1) x (2) w w x (2) x (1) x (1) x ( ) x (1) x ( ) Fig. 2 . Ensemble of 3-point quadrature rules on
Γ = [0 , found via designed quadrature( d = r = 2 ). Each three-point nodal configuration has nodes connected with blue lines,forming a triangle. Left: initial guesses provided to the algorithm. Right: converged designedquadrature rules. Bottom: Nodal configurations on Γ . Top: Weight values plotted as z -coordinates. in [24], we can explicitly compute this as L (Λ T r ) = (cid:12)(cid:12) Λ T (cid:98) r/ (cid:99) (cid:12)(cid:12) = (cid:18) d + (cid:98) r/ (cid:99) d (cid:19) (33)Independently, we computed designed quadratures for r = 2 and r = 3 to confirmthat the number of nodes for different dimensions d coincides with d + 1 and 2 d ,respectively, as determined in [41] (not shown). Also, for r = 5 and d = 3 , d = 4 for different orders r and we setthe tolerance to (cid:15) = 10 − in this example. It should be noted that these quantitativemetrics can vary depending on the random initialization and regularization parametersthroughout the algorithm however provide a useful holistic measure for the method’sperformance.To illustrate how the regularization parameter λ is chosen, we show the singularvalues and regularization parameter choice for the case r = 5 , d = 7. Figure 4 showsthe regularization parameter selection for an iteration in the middle of the procedure.The regularized parameter is selected as λ = 10 by investigating the spectrum ofsingular values and its L-curve. This manuscript is for review purposes only.
ESIGNED QUADRATURE Table 1
Number of nodes n sparse grids and designed quadratures on total degree spaces Λ T r on Γ =[0 , d . Top: fixed d = 3 for various r . Bottom: Fixed r = 5 for various d . The results in thetop-half of this table can be compared with Table 4 in [53]. Our quadrature rules have smaller orequal size compared with the results in [53], with the exception of r = 8 where we report a 43-pointrule instead of a 42-point rule in [53]. d = 3 , r r = 5 , d Table 2
Performance of the scheme with respect to number of nodes and iterations, CPU time andresidual norm for d = 4 and various orders r with total order index set. d = 4 , r (cid:12)(cid:12) Λ T (cid:98) r/ (cid:99) (cid:12)(cid:12) || ˜ R || In practice, one could fix λ as a function of (cid:107) R (cid:107) (or (cid:107) (cid:101) R (cid:107) ). In Figure 3, we have λ = 10 with (cid:107) R (cid:107) = 40. Then, for example, one could take λ = 50 for 200 ≤ || R || ≤
500 and λ = 10 for 20 ≤ || R || ≤ a priori tabulation could be fixed fora variety of ( d, r ) values. Designed quadrature rulescan be used to construct polynomial interpolants. Suppose we have a designed quadra-ture rule ( X , w ) of size n that matches moments for indices on Λ (up to the tolerance (cid:15) ), and assume n = | Λ | . For continuous function f , let I ( f ) denote the uniqueinterpolant of f from Π Λ at the locations X . Lebesgue’s lemma states || f − I ( f ) || ∞ ≤ ( L + 1) inf p ∈ Π Λ || f − p || ∞ , L = sup (cid:107) h (cid:107) ∞ =1 (cid:107)I ( h ) (cid:107) ∞ , where (cid:107)·(cid:107) ∞ is the maximum norm on Γ, and the supremum is taken over all functions h continuous on Γ. The constant L is the Lebesgue constant; small values indicatethat interpolants are comparable to the best approximation measured in the maximumnorm [27]. The Lebesgue constant can be computed explicitly: The interpolant I ( f )can be expressed as I ( f )( x ) = n (cid:88) j =1 (cid:96) j ( x ) f ( x j ) , (cid:96) j ( x k ) = δ j,k , Designed quadrature rules achieve n < | Λ | , but in this section we will enforce n = | Λ | for thepurposes of forming an interpolant. This manuscript is for review purposes only. V. KESHAVARZZADEH, R. M. KIRBY, AND A. NARAYAN r n Lower boundDes. Quad.SG-KPProd. ruleSG-GQ d n Lower boundDes. Quad.SG-KPSG-GQProd. rule
Fig. 3 . Number of nodes for fixed d = 3 (left), and fixed r = 5 (right) for total or-der index set Λ T r . The lower bound is given in (33) , “SG-KP” and “SG-GQ” are sparsegrid constructions using nested Kronrod-Patterson and non-nested Gauss quadrature rules,respectively. where the (cid:96) j are the cardinal interpolation functions. The Lebesgue function L n andthe Lebesgue constant L are, respectively, L n ( x ) = n (cid:88) j =1 | (cid:96) j ( x ) | , L = (cid:107) L n (cid:107) ∞ Finding a set of points with minimal Lebesgue constant is not trivial. In d = 2 di-mensions the Padua points are essentially the only explicitly constructible set of nodeswith provably minimal growth of Lebesgue constant on total degree spaces [5]. Tocompare designed quadrature with Padua points, we consider degree-5 Padua points,yielding dim Π Λ T = 21. These points along with associated quadrature weights inte-grate polynomials in Π Λ T exactly with respect to the product Chebyshev weight ω [5]. With designed quadrature we are able to find 17 <
21 nodes and weights thatintegrate polynomials in Π Λ T exactly. However, for the purposes of interpolationin this section, we enforce n = 21 nodes in the designed quadrature framework. Toinitialize the design we start from nodes that are close to Padua nodes. The Lebesguefunction L n ( x ) for both cases are shown in Figure 5. The Lebesgue constant for Paduapoints and designed quadrature are L = 4 . L = 5 . L suggest that the designed quadrature points and the Paduapoints are of comparable quality in terms of constructing interpolants. However, wereiterate that for quadrature we can use fewer nodes (17) than the Padua points (21). The formulation of designed quadrature allowsΓ and ω to be of relatively general form, but we can construct quadrature rules ineven more exotic situations. Let Γ = [ − , with ω the uniform weight. Instead ofenforcing x j ∈ Γ in (21), we enforce x j ∈ (cid:101) Γ, where (cid:101) Γ ⊂ Γ is a “U” shape, mimickingthe logo of the University of Utah; see Figure 6, left.The penalty function for this problem has the same quadratic form as thosediscussed in Section 3.1 and separate penalties are considered for violations in both x (1) and x (2) directions. For example, we can model the infeasible rectangular region S between the two ascenders of the U with non-zero penalty in x (1) direction and This manuscript is for review purposes only.
ESIGNED QUADRATURE i -0.8-0.6-0.4-0.20 ∆ l og ( σ i ) i -3 -2 -1 σ i
10 20 40 50 70 || e J ∆ d λ − e R || -5 || ∆ d λ || k -10 -5 M agn i t ude || e R || η Fig. 4 . Singular values of (cid:101) J with the chosen regularization parameter (top left), finitedifference on log of singular values and the chosen singular value index (top right), L -curvefor the given (cid:101) J and (cid:101) R : the point on the L -curve corresponding to the selected regularizationparameter λ is indicated with the red circle (bottom left), convergence of the scheme for d = 5 , r = 7 (bottom right). zero penalty in x (2) direction as(34) (cid:40) S : (0 ≤ | x (1)1 | ≤ . ∩ ( − . ≤ x (2) ≤ . P = ( x (1) − . , P = 0A similar method can be used to penalize the semicircular region below the rectanglewhere violations in both directions are penalized. The total penalty for infeasibleregions then involves both P and P e.g. P = 10 √ P + P which is shown in Figure6, right. We compute a designed quadrature rule for total degree r = 2, achievingresidual tolerance of (cid:15) = 0 . n = 150. We need a relatively large numberof nodes, and achieve only a relatively large tolerance (compared to (cid:15) = 10 − forprevious examples). This is due to the difficulty of this problem: we want nodesto lie in (cid:101) Γ but want to achieve integration over Γ. We expect that convergence forlarger r will require many iterations and may not be able to achieve arbitrarily smalltolerances. To demonstrate the capability of de-signed quadrature for integration in high dimensions, we consider d = 100 with hy- This manuscript is for review purposes only. V. KESHAVARZZADEH, R. M. KIRBY, AND A. NARAYAN
Fig. 5 . Contour plots of Lebesgue function L n for Padua points (left) and designedquadrature (right) -1 -0.5 0 0.5 1 x (1) -1-0.500.51 x ( ) x (2)
01 0.5 x (1) -0.50 -0.5 -1-1 -1.5-1.520 P Fig. 6 . Designed Quadrature for uniform weight and d = r = 2 with “U” shape indicatingthe University of Utah (left), the penalty function used in the scheme (right). perbolic cross index set Λ H r . Figure 7 (top) shows different slices of the d = 100nodal configuration generated by designed quadrature for the uniform weight onΓ = [ − , and r = 4, for which we have n = 106 and | Λ H | = 5351.Figure 7 (bottom) shows the behavior of designed quadrature weights with respectto the Euclidean norm of the nodes (distance to the origin) for ω the Gaussian weighton Γ = R for total order r = 2 and hyperbolic cross orders r = 3 ,
4. As expected theweights decay as the node norms increase. We find n = 101 for all these quadratures,again confirming the optimal n = d + 1 size for total order r = 2 [54]. It is alsointeresting to note that the minimum Euclidean norm of nodes for these cases aresomewhat equal viz || x || = 8 . , . , .
85 respectively.Computing designed quadratures in high dimensions reveals computational chal-lenges that are not present in small-to-moderate dimensions: since the nonlinearsystem is quite large, we do not perform the SVD of Jacobian in each iteration. In-stead, we regularize the pseudo-inverse matrix directly and compute the Newton stepas ∆ d = ( J T J + λ I ) − J T R . The parameter λ can be selected based on the residual This manuscript is for review purposes only.
ESIGNED QUADRATURE d = 100 took ∼
100 iterations and less than 30 minutes on a personaldesktop in MATLAB. -1 0 1 x (1) -1-0.500.51 x ( ) -1 0 1 x (10) -1-0.500.51 x ( ) -1 0 1 x (50) -1-0.500.51 x ( ) ||x i || w i ||x i || w i ||x i || w i Fig. 7 . Different two dimensional slices of d = 100 designed quadrature for uniformweight and hyperbolic cross order r = 4 (top), formation of weights with respect to Euclideannorm of nodes for Gaussian weight in d = 100 (bottom): total order r = 2 (left) hyperboliccross order r = 3 (middle) and r = 4 (right). To inves-tigate the performance of the high dimensional quadrature points we compute themean and variance of compliance indicative of the elastic energy for a solid cantileverbeam with uncertain material properties.The compliance for the spatial domain Ω reads C = (cid:90) Ω f ud Ωwhere u is the displacement and f is the surface load on the structure. To find u ,the equation of motion in linear elasticity ∇ . σ + f = 0 where σ is the stress tensorand ∇ is the divergence operator is solved via Finite Element Method. The globaldisplacement is characterized with n e finite elements u = n e (cid:88) i =1 u i Ψ i where Ψ i are finite element shape functions and u i are nodal displacements. Thenodal displacements U = { u i } n e i =1 are solution of a linear system KU = F (stems This manuscript is for review purposes only. V. KESHAVARZZADEH, R. M. KIRBY, AND A. NARAYAN from the equation of motion) where K = (cid:90) Ω ∂ Ψ T ∂ x C ∂ Ψ ∂ x d Ω , F = (cid:90) Ω f Ψ d Ωwith C being an elasticity matrix. We consider the plane stress condition in thisexample, hence for our two dimensional problem C = E − ν ν ν − ν where E is the modulus of elasticity and ν is the Poisson’s ratio.The beam geometry is shown in Figure 8, and is modeled with 100 standardsquare finite elements, where each element has lognormal modulus of elasticity as E i = 10 − +exp( ξ ( i ) ). The random variables ξ ( i ) | i =1 are independent standard randomnormal variables and the Poisson’s ratio is ν = 0 .
3. Our goal is to compute first- andsecond-order statistics of the compliance; these statistics are integrals with respect tothe 100 variables ξ ( i ) , and so we approximate these statistics via designed quadrature.For comparison against designed quadrature we use quasi-Monte Carlo (QMC)samples of size n = 64 , , ,
512 and n = 1024, where we treat the latter asthe exact solution. The QMC samples are generated on [0 , , and are mapped to R using inverse transform sampling for a 100-dimensional standard normal randomvector. We have chosen the number of QMC samples so that they almost match thenumber of designed quadrature nodes computed from the index sets (i) total orderwith r = 2 and n = 101 nodes, (ii) the index set Λ H ∪ Λ with n = 155 nodes,where Λ contains pairwise interactions of maximum univariate order 2, and (iii) theindex set Λ H ∪ Λ with n = 255 nodes, where Λ contains pairwise interactions ofmaximum univariate order 3.Figure 8 compares errors in the computed mean and standard deviation of thecompliance for QMC versus designed quadrature; we see that designed quadratureachieves significantly better errors. We believe QMC would be more effective if theproblem involved many more variables, larger index sets, and/or non-smooth quanti-ties of interest, as in those cases the number of designed quadrature is prohibitivelylarge and finding a suitable quadrature rule is computationally challenging. Our final example utilizes polynomial chaos (PC) methods [17] to build surrogates fortopology optimization under geometric uncertainty [25]. Figure 9 shows the flowchartfor design optimization under uncertainty. To build PC surrogates at each designiteration, n Finite Element Analysis (FEA) and sensitivity analyses are performed inorder to quantify the uncertainty associated with random variables ξ ( i ) as in the lastexample. This is the most costly step in the design process, and hence small n canresult in significant computational savings.The perturbation in the boundary of topology interfaces Z are modeled via a This manuscript is for review purposes only.
ESIGNED QUADRATURE
64 128 256 512 n -4 -3 -2 -1 R e l a t i v e e rr o r QMC-MeanQMC-StdDQ-MeanDQ-Std
Fig. 8 . Finite element discretization for a linear elastic cantilever beam with randomelastic modulus (left) convergence of compliance mean and variance with quasi-Monte Carlo(QMC) samples and designed quadrature (DQ) in d = 100 variables (right). Design Optimization under UncertaintyGradient-Based Optimizer θ θ ∗ Forward Model (FEA)& Sensitivity Analysis f ( θ , ξ i ) | n i =1 ∂f ( θ , ξ i ) ∂ θ | n i =1 Surrogate Model (PCE) ⇓ Performance Metrics& Their Sensitivities µ ( θ ) , σ ( θ ) ∂µ ( θ ) ∂ θ , ∂σ ( θ ) ∂ θ RandomVariables ξ i | n i =1 DesignVariables θ Fig. 9 . Design optimization under uncertainty flowchart
Karhunen-Lo`eve random field with d = 4 significant modes as(35) Z ( x, ξ ) = (cid:88) i =1 (cid:112) λ i γ ( i ) ( x ) ξ ( i ) . where ξ ( i ) ∈ U [ −√ , √
3] are independent uniform random variables. Sparse gridsbuilt from nested rules were utilized in [25] to develop a surrogate for total degree r = 3in d = 4 dimensions. This requires a quadrature rule that can accurately integratepolynomials up to order r = 6 (see Proposition 2.3). A standard construction of sparsegrid rules yields odd orders of polynomial accuracy and hence a quadrature rule for r = 7 with n = 81 nodes is used. We observe that 33 out of these 81 nodes havenegative weights [22]. On the other hand we use designed quadrature constrained tointegrate polynomials up to degree r = 6 and compute n = 43 nodes, almost half( ∼ n = 641 points ( r = 13)as the “true” solution. The mean and standard deviation for the true solution This manuscript is for review purposes only. V. KESHAVARZZADEH, R. M. KIRBY, AND A. NARAYAN
118 119 120 121 122
Compliance P D F Sp. GridDes. Quad.
Fig. 10 . Robust topology design of the MBB beam (left) Compliance PDF with SparseGrid and designed quadrature (right) are µ true = 120 . , σ true = 0 . e µ = | ( µ − µ true ) /µ true | and in standard deviation e σ = | ( σ − σ true ) /σ true | are listed in Table 3. We achieve higher accuracy withdesigned quadrature at nearly half the cost. Figure 10 also visually compares theprobability density function (PDF) of compliance for both cases, and no substantialdifference is observed. Table 3
Mean and standard deviation estimation for the robust topology design.
Quadrature Rule µ e µ σ e σ CostSparse Grid 120.1326 2.44e-04 0.7836 2.16e-03 81 SimulationsDesigned Quadrature 120.1028 3.33e-06 0.7854 1.27e-04 43 Simulations
5. Concluding Remarks.
We present a systematic approach, designed quadra-ture, for computing multivariate quadrature rules in generic settings. The frameworkuses penalty methods in constrained optimization to ensure positivity of the weightsand feasible locations for the nodes. The Gauss-Newton algorithm is used to performminimization of the penalty-augmented objective function. L regularization is uti-lized to treat ill-conditioned systems encountered during Newton step updates. Onregular domains such as hypercubes, our designed quadrature results in considerablyfewer nodes (and guaranteed positive weights) compared to alternative multivariatequadrature rules, such as sparse grids, and hence is promising for computational sci-ence and engineering involving expensive simulations. When applied to a benchmarkrobust topology optimization problem, designed quadrature reduces requisite cost bynearly half compared with sparse grid rules, and achieves higher accuracy. REFERENCES[1]
N. Agarwal and N. R. Aluru , A domain adaptive stochastic collocation approach foranalysis of MEMS under uncertainties
I. Babuka, F. Nobile, and R. Tempone , A Stochastic Collocation Method for Elliptic Par-tial Differential Equations with Random Input Data , SIAM Review, 52 (2010), pp. 317–355, https://doi.org/10.1137/100786356, http://epubs.siam.org.libproxy.umassd.edu/doi/abs/10.1137/100786356 (accessed 2014-06-23).
This manuscript is for review purposes only.
ESIGNED QUADRATURE Table 4
Designed nodes and weights for uniform wight function associated with d = 4 , r = 6 . x (1) x (2) x (3) x (4) w [3] C. Bernardi and Y. Maday , Spectral methods
D. Bertsekas , Nonlinear programming , Athena Scientific, Second Edition, (2008).[5]
L. Bos, M. Caliari, M. Vianello, S. De Marchi, and Y. Xu , Bivariate lagrange interpolationat the padua points: the generating curve approach , Journal of Approximation Theory, 143(2006), pp. 15–25.[6]
S. Boyd and L. Vandenberghe , Convex optimization , Cambridge University Press, (2004).[7]
H. Bungartz and M. Griebel , Sparse grids , Acta Numerica, 13 (2004), pp. 147 – 269.[8]
H.-J. Bungartz and S. Dirnstorfer , Multivariate Quadrature on Adaptive Sparse Grids
M. Caliari, S. De Marchi, and M. Vianello , Bivariate polynomial interpolation on thesquare at new nodal sets , Applied Mathematics and Computation, 165 (2005), pp. 261–274, https://doi.org/10.1016/j.amc.2004.07.001.[10]
C. Canuto, M. Y. Hussaini, A. Quarteroni, and T. A. Zang , Spectral Methods: Funda-mentals in Single Domains , Springer, Berlin ; New York, 1st ed. 2006. corr. 4th printing2010 edition ed., Sept. 2011.[11]
C. Canuto, M. Y. Hussaini, A. Quarteroni, and T. A. Zang , Spectral Methods: Evolu-tion to Complex Geometries and Applications to Fluid Dynamics , Springer, Berlin, 2007edition ed., Nov. 2014.[12]
A. Cohen, R. DeVore, and C. Schwab , Convergence Rates of Best N-term Galerkin Ap-proximations for a Class of Elliptic sPDEs , Foundations of Computational Mathematics,10 (2010), pp. 615–646, https://doi.org/10.1007/s10208-010-9072-2, https://link.springer.com/article/10.1007/s10208-010-9072-2 (accessed 2017-10-18).[13]
P. G. Constantine, M. S. Eldred, and E. T. Phipps , Sparse pseudospectral approximationmethod , Computer Methods in Applied Mechanics and Engineering, 229 (2012), pp. 1–12,https://doi.org/10.1016/j.cma.2012.03.019.
This manuscript is for review purposes only. V. KESHAVARZZADEH, R. M. KIRBY, AND A. NARAYAN[14]
P. Davis and P. Rabinowitz , Methods of numerical integration , Courier Corporation, 2 (2007).[15]
W. Gautschi , Construction of Gauss-Christoffel quadrature formulas , Mathematics of Com-putation, 22 (1968), pp. 251–270.[16]
T. Gerstner and M. Griebel , Numerical integration using sparse grids , Numerical Algo-rithms, 18 (1998), pp. 209–232, https://doi.org/10.1023/A:1019129717644.[17]
R. Ghanem and P. Spanos , Stochastic finite elements: A spectral approach , Dover publica-tions, (2002).[18]
G. Golub and J. Welsch , Calculation of gauss quadrature rules , Mathematics of Computa-tion, 23 (1969), pp. 221 – 230.[19]
G. H. Golub and C. F. V. Loan , Matrix Computations Johns Hopkins Studies in Mathemat-ical Sciences , The Johns Hopkins University Press, 3rd ed., Oct. 1996.[20]
P. Hansen , Rank-deficient and discrete ill-posed problems , SIAM, Philadelphia, (1998).[21]
P. Hansen and D. OLeary , The use of the L-curve in the regularization of discrete ill-posedproblems , SIAM Journal on Scientific Computing, 14 (1993), pp. 1487–1503.[22]
F. Heiss and V. Winschel , Quadrature on sparse grids
F. Heiss and V. Winschel , Likelihood approximation by numerical integration on sparse grids ,Journal of Econometrics, 144 (2008), pp. 62 – 80.[24]
J. Jakeman and A. Narayan , Generation and application of multivariate polynomial quadra-ture rules , arXiv:1711.00506 [math.NA], (2017).[25]
V. Keshavarzzadeh, F. Fernandez, and D. Tortorelli , Topology optimization under un-certainty via non-intrusive polynomial chaos expansion , Computer Methods in AppliedMechanics and Engineering, 318 (2017), pp. 120–147.[26]
M. Liu, Z. Gao, and J. S. Hesthaven , Adaptive sparse grid algorithms with applications toelectromagnetic scattering under uncertainty , Applied Numerical Mathematics, 61 (2011),pp. 24–37, https://doi.org/10.1016/j.apnum.2010.08.002.[27]
D. Lubinsky , A survey of weighted polynomial approximation with exponential weights , Surveysin Approximation Theory, 3 (2007), pp. 1–105.[28]
D. Luenberger and Y. Ye , Linear and nonlinear programming , International Series in Oper-ations Research and Management Science, (2008).[29]
J. Ma, V. Rokhlin, and S. Wandzura , Generalized gaussian quadrature rules for systems ofarbitrary functions , SIAM Journal on Numerical Analysis, 33 (1996), pp. 971–996.[30]
S. E. Mousavi, H. Xiao, and N. Sukumar , Generalized Gaussian quadrature rules on arbitrarypolygons , International Journal for Numerical Methods in Engineering, 82 (2010), pp. 99–113, https://doi.org/10.1002/nme.2759.[31]
A. Narayan and J. Jakeman , Adaptive Leja Sparse Grid Constructions for StochasticCollocation and High-Dimensional Approximation , SIAM Journal on Scientific Comput-ing, 36 (2014), pp. A2952–A2983, https://doi.org/10.1137/140966368. arXiv:1404.5663[math.NA].[32]
A. Narayan, J. Jakeman, and T. Zhou , A christoffel function weighted least squares algorithmfor collocation approximations , Mathematics of Computation, 86 (2017), pp. 1913–1947.[33]
H. Niederreiter , Random number generation and quasi-monte carlo methods , Society forIndustrial and Applied Mathematics, (1992).[34]
A. B. Owen , Quasi-monte carlo sampling , Monte Carlo Ray Tracing: SIGGRAPH 2003 Course44, (2003), pp. 69–88.[35]
T. Patterson , The optimum addition of points to quadrature formulae , Mathematics of Com-putation, 22 (1968), pp. 847–856.[36]
E. Ryu and S. Boyd , Extensions of gauss quadrature via linear programming , Foundations ofComputational Mathematics, 15 (2015), pp. 953–971.[37]
I. H. Sloan and H. Wozniakowski , When Are Quasi-Monte Carlo Algorithms Efficient forHigh Dimensional Integrals? , Journal of Complexity, 14 (1998), pp. 1–33, https://doi.org/10.1006/jcom.1997.0463.[38]
S. Smolyak , Quadrature and interpolation formulas for tensor products of certain classes offunctions , Soviet Mathematics Doklady, 4 (1963), pp. 240–243.[39]
J. Stoer and R. Bulirsch , Introduction to numerical analysis , Springer-Verlag New York, 12(2002).[40]
A. Stroud , Some fifth degree integration formulas for symmetric regions II , Numerische Math-ematik, 9 (1967), pp. 460–468.[41]
A. Stroud , Approximate calculation of multiple integrals , Englewood Cliffs, N.J., Prentice-Hall, (1971).[42]
A. H. Stroud , Remarks on the Disposition of Points in Numerical Integration Formulas ,Mathematical Tables and Other Aids to Computation, 11 (1957), pp. 257–261.[43]
A. H. Stroud , Numerical Integration Formulas of Degree Two , Mathematics of Computation,
This manuscript is for review purposes only.
ESIGNED QUADRATURE
14 (1960), pp. 21–26, https://doi.org/10.2307/2002981.[44]
G. Szeg¨o , Orthogonal Polynomials , American Mathematical Soc., 4th ed., 1975.[45]
M. A. Taylor, B. A. Wingate, and L. P. Bos , A Cardinal Function Algorithm for ComputingMultivariate Quadrature Points , SIAM Journal on Numerical Analysis, 45 (2007), pp. 193–205, https://doi.org/10.1137/050625801.[46]
M. A. Taylor, B. A. Wingate, and R. E. Vincent , An Algorithm for Computing FeketePoints in the Triangle , SIAM Journal on Numerical Analysis, 38 (2000), pp. 1707–1720,https://doi.org/10.1137/S0036142998337247.[47]
R. Tibshirani , Regression Shrinkage and Selection via the Lasso , Journal of the Royal Statis-tical Society. Series B (Methodological), 58 (1996), pp. 267–288.[48]
M. van Barel, M. Humet, and L. Sorber , Approximating optimal point configurations formultivariate polynomial interpolation , Electronic Transactions on Numerical Analysis, 42(2014), pp. 41–63.[49]
E. Van den berg and M. Friedlander , Spgl1: A solver for large-scale sparse reconstruction
E. Van den berg and M. Friedlander , Probing the Pareto frontier for basis pursuit solutions ,SIAM Journal on Scientific Computing, 31 (2008), pp. 890–912.[51]
E. Van den berg and M. Friedlander , Sparse optimization with least-squares constraints ,SIAM Journal on Optimization, 21 (2011), pp. 1201–1229.[52]
G. Wasilkowski and H. Wozniakowski , Explicit cost bounds of algorithms for multivariatetensor product problems , Journal of Complexity, 11 (1995), pp. 1 – 56.[53]
H. Xiao and Z. Gimbutasb , A numerical algorithm for the construction of efficient quadraturerules in two and higher dimensions , Computers and Mathematics with Applications, 59(2010), pp. 663–676.[54]
D. Xiu , Numerical integration formulas of degree two , Applied Numerical Mathematics, 58(2008), pp. 1515–1520, https://doi.org/16/j.apnum.2007.09.004.[55]
D. Xiu and J. S. Hesthaven , High-Order Collocation Methods for Differential Equationswith Random Inputs , SIAM Journal on Scientific Computing, 27 (2005), pp. 1118–1139,https://doi.org/10.1137/040615201.[56]
Y. Xu , A characterization of positive quadrature formulae , Mathematics of Computation, 62(1994), pp. 703–718, https://doi.org/10.1090/S0025-5718-1994-1223234-0., Mathematics of Computation, 62(1994), pp. 703–718, https://doi.org/10.1090/S0025-5718-1994-1223234-0.