AN EFFICIENT AND FAST SPARSE GRID ALGORITHM FOR HIGH-DIMENSIONAL NUMERICAL INTEGRATION HUICONG ZHONGAND XIAOBING FENG

2025-04-30 0 0 872.83KB 28 页 10玖币

侵权投诉

AN EFFICIENT AND FAST SPARSE GRID ALGORITHM FOR

HIGH-DIMENSIONAL NUMERICAL INTEGRATION

HUICONG ZHONG∗AND XIAOBING FENG†

Abstract. This paper is concerned with developing an eﬃcient numerical algorithm for fast

implementation of the sparse grid method for computing the d-dimensional integral of a given func-

tion. The new algorithm, called the MDI-SG (multilevel dimension iteration sparse grid) method,

implements the sparse grid method based on a dimension iteration/reduction procedure, it does

not need to store the integration points, neither does it compute the function values independently

at each integration point, instead, it re-uses the computation for function evaluations as much as

possible by performing the function evaluations at all integration points in cluster and iteratively

along coordinate directions. It is showed numerically that the computational complexity (in terms

of CPU time) of the proposed MDI-SG method is of polynomial order O(N d3) or better, compared

to the exponential order O(N(log N)d−1) for the standard sparse grid method, where Ndenotes

the maximum number of integration points in each coordinate direction. As a result, the proposed

MDI-SG method eﬀectively circumvents the curse of dimensionality suﬀered by the standard sparse

grid method for high-dimensional numerical integration.

Key words. Sparse grid(SG), multilevel dimension iteration (MDI), high-dimensional integra-

tion, numerical quadrature rules, curse of dimensionality.

AMS subject classiﬁcations. 65D30, 65D40, 65C05, 65N99

1. Introduction. With rapid developments in nontraditional applied sciences

such as mathematical ﬁnance [11], image processing [15], economics [1], and data

science [24], there is an ever increasing demand for eﬃcient numerical methods for

computing high-dimensional integration which also becomes crucial for solving some

challenging problems. Numerical methods (or quadrature rules) mostly stem from

approximating the Riemann sum in the deﬁnition of integrals, hence, they are grid-

based. The simplest and most natural approach for constructing numerical quadrature

rules in high dimensions is to apply the same 1-d rule in each coordinate direction,

this then leads to tensor-product (TP) quadrature rules. It is well known (and easy

to check) that the number of integration points (and function evaluations) grows

exponentially in the dimension d, such a phenomenon is known as the curse of di-

mensionality (CoD). Mitigating or circumventing the CoD has been the primary goal

when it comes to constructing eﬃcient high-dimensional numerical quadrature rules.

A lot of progress has been made in this direction in the past ﬁfty years, this includes

sparse grid (SG) methods [4,10,11,7], Monte Carlo (MC) methods [5,21], Quasi-

Monte Carlo (QMC) methods [6,13,14,16,27], deep neural network (DNN) methods

[8,12,17,25,28]. To some certain extent, those methods are eﬀective for computing

integrals in low and medium dimensions (i.e., d≲100), but it is still a challenge for

them to compute integrals in very high dimensions (i.e., d≈1000).

This is the second installment in a sequel [9] which aims at developing fast nu-

merical algorithms for high-dimensional numerical integration. As mentioned above,

the straightforward implementation of the TP method will evidently run into the

CoD dilemma. To circumvent the diﬃculty, we proposed in [9] a multilevel dimension

iteration algorithm (called MDI-TP) for accelerating the TP method. The ideas of

the MDI-TP algorithm are to reuse the computation of function evaluation as much

∗School of Mathematics and Statistics, Northwestern Polytechnical University, Xi’an, Shaanxi

710129, China (huicongzhong@mail.nwpu.edu.cn).

†Department of Mathematics, The University of Tennessee, Knoxville, TN, 37996

(xfeng@utk.edu).

2HUICONG ZHONG AND XIAOBING FENG

as possible in the tensor product method by clustering computations, which allows

an eﬃcient and fast function evaluations at integration points together, and to do

clustering by a simple dimension iteration/reduction strategy, which is possible be-

cause of the lattice structure of the TP integration points. Since the idea of the MDI

strategy essentially applied to any numerical integration rule whose integration points

have a lattice-like structure, this indeed motivates the work of this paper by applying

the MDI idea to accelerate the sparse grid method.

The sparse grid (SG) method, which was ﬁrst proposed by Smolyak in [26], only

uses a (small) subset of the TP integration points while still maintains a comparable

accuracy of the TP method. As mentioned earlier, the SG method was one of few suc-

cessful numerical methods which can mitigate the CoD in high-dimensional computa-

tion, including computing high-dimensional integration and solving high-dimensional

PDEs [4,10,11,7]. The basic idea in application to high-dimensional numerical inte-

gration stems from Smolyak’s general method for multivariate extensions of univariate

operators. Based on this construction, the midpoint rule [2], the rectangle rule [22],

the trapezoidal rule [3], the Clenshaw-Curtis rule [19,20], and the Gaussian-Legendre

rule [18,23] have been used as a one-dimensional numerical integration method. A

multivariate quadrature rule is then constructed by forming the TP method of each of

these one-dimensional rules on the underlying sparse grid. Like TP method, the SG

method is quite general and easy to implement. But unlike the TP method, its com-

putational cost is much lower because the number of its required function evaluations

grows exponentially with much smaller base.

The goal of this paper is to apply the MDI strategy to the SG method, the

resulting algorithm, which is called the SG-MDI algorithm, provides a fast algorithm

for an eﬃcient implementation of the SG method. The SG-MDI method incorporates

the MDI nested iteration idea into the sparse grid method which allows the reuse of the

computation in the function evaluations at the integration points as much as possible.

This saving signiﬁcantly reduces the overall computational cost for implementing the

sparse grid method from an exponential growth to a low-order polynomial growth.

The rest of this paper is organized as follows. In section 2, we ﬁrst brieﬂy recall

the formulation of the sparse grid method and then introduce our SG-MDI algorithm

in two and three dimensions to explain the main ideas of the algorithm as well as

generalize it to arbitrary dimensions. In section 4, we present various numerical

experiments to test the performance of the proposed SG-MDI algorithm and compare

its performances to the standard SG method and the classical MC method. These

numerical tests show that the SG-MDI algorithm is much faster in low and medium

dimensions (i.e. d≲100) and in very high dimensions (i.e. d≈1000). It still works

even when the MC method fails to compute. In section 5, we provide numerical

tests to measure the inﬂuence of parameters in the proposed SG-MDI algorithm,

including dependencies on the choice of underlying 1-d quadrature rule, the precision

layer, and the choice of iteration step size. In section 6we numerically estimate the

computational complexity of the SG-MDI algorithm. This is done by using standard

regression technique to discover the relationship between CPU time and dimension.

It is showed that the CPU time grows at most in a polynomial order O(d3N), where d

and Nstand for respectively the dimension of integration domain and the maximum

number of integration points in each coordinate direction. As a result, the SG-MDI

algorithm can easily compute intermediate and very high-dimensional integrals on

common desktop computers. Finally, the paper is concluded with some concluding

remarks given in section 7.

A FAST SG-MDI ALGORITHM FOR HIGH-DIMENSIONAL INTEGRATIPON 3

2. Preliminaries. In this paper, f(x) : ¯

Ω→Rdenote a generic continuous

function on ¯

Ω⊂Rdfor d >> 1, then f(x) has pointwise values at every x=

(x1, x2,··· , xd)∈¯

Ω.Without loss of the generality, we set Ω := [−1,1]dand con-

sider the basic and fundamental problem of computing the integral

(2.1) Id(f) := ZΩ

f(x)dx.

As mentioned in Section 1, the goal of this paper is to develop a fast algorithm for

computing above integral based on the sparse grid methodology. To that end, below

we brieﬂy recall the necessary elements of sparse grid methods.

2.1. The sparse grid method. We now recall the formulation of the sparse

grid method for approximating (2.1) and its tensor product reformulation formula

which will be crucially used later in the formulation of our fast SG-MDI algorithm.

For each positive integer index l≥1, let nlbe a positive integer which denotes

the number of grid points at level land

(2.2) Γl

1:= xl

1< xl

2<··· < xl

nl⊂[−1,1]

denote a sequence of the level lgrid points in [−1,1]. The grid sets {Γl

1}is said to be

nested provided that Γl

1⊂Γl+1

1. The best known example of the nested grids is the

following dyadic grids:

(2.3) Γl

1=nxl

i:= i

2l−1:i=−2l−1,··· ,0,1,··· ,2l−1o.

For a given positive integer q, the tensor product Gq

d:= Γq

1×Γq

1× ··· × Γq

then yields a standard tensor product grid on Ω = [−1,1]d. Notice that here the

qth level is used in each coordinate direction. To reduce the number of grid points

in Gq

d, the sparse grid idea is to restrict the total level to be qin the sense that

q=l1+l2+··· +ld, where liis the level used in the ith coordinate direction, its

corresponding tensor product grid is Γl1

1×Γl2

1×···×Γld

1. Obviously, the decomposition

q=l1+l2+···+ldis not unique, so all such decomposition must be considered. The

union

(2.4) Γq

d:= [

l1+···+ld=q

Γl1

1× ··· × Γld

then yields the famous Smolyak sparse grid on Ω = [−1,1]d(cf. [7]). We remark

that the underlying idea of going from Gq

dto Γq

dis exactly same as going from Qq(Ω)

to Pq(Ω), where Qqdenotes the set of polynomials whose degrees in all coordinate

directions not exceeding qand Pqdenotes the set of polynomials whose total degrees

not exceeding q.

After having introduced the concept of sparse grids, we then can deﬁne the sparse

grid quadrature rule. For a univariate function gon [−1,1], we consider done-

dimensional quadrature formula

(2.5) Jli

1(g) :=

nli

j=1

wli

jg(xli

j), i = 1,··· , d,

where {xli

j, j = 1,··· , nli}and {wli

j, j = 1,··· , nli}denote respectively the integra-

tion points/nodes and weights of the quadrature rule, and nlidenotes the number of

4HUICONG ZHONG AND XIAOBING FENG

integration points in the ith coordinate direction in [−1,1]. Deﬁne

l:= [l1, l2,··· , ld],|l|:= l1+l2+··· +ld,

Nd,s := nl= [l1,··· , ld] : |l|=s, s ≥do.

For example N2,4=[1,3],[2,2],[3,1].

Then, the sparse grid quadrature rule with accuracy level q∈Nfor d-dimensional

integration (2.1) on [−1,1]dis deﬁned as (cf. [11])

(2.6) Qq

d(f) := X

q−d+1≤|l|≤q

(−1)q−|l|d−1

q− |l|X

l∈Nd,|l|Jl1

1⊗ Jl2

1⊗ ··· ⊗ Jld

1f,

where

(2.7) Jl1

1⊗ Jl2

1⊗ ··· ⊗ Jld

1f:=

nl1

j1=1 ···

nld

jd=1

wl1

j1···wld

jdfxl1

j1,··· , xld

jd.

We note that each term (Jl1

1⊗ Jl2

1⊗ ··· ⊗ Jld

1)fin (2.7) is the tensor product

quadrature rule which uses nliintegration points in ith coordinate direction. To write

d(f) more compactly, we set

d:= X

q−d+1≤|l|≤q

nl1. . . nld,

which denotes the total number of integration points in Ω = [−1,1]d. Let wk, k =

1,··· , nq

ddenote the corresponding weights and deﬁne the bijective mapping

xk:k= 1,··· , nq

d−→ (xl1

j1,··· , xld

jd) : ji= 1,··· , nli, q −d+ 1 ≤ |l| ≤ q.

Then, the sparse grid quadrature rule Qq

d(f) can be rewritten as

(2.8) Qq

d(f) =

k=1

wkf(xk).

We also note that some weights wkmay become negative even though the one-

dimensional weights wl1

j1,··· , wld

jdare positive. Therefore, it is no longer possible

to interpret Qq

d(f) as a discrete probability measure. Moreover, the existence of

negative weights in (2.8) may cause numerical cancellation, hence, loss of signiﬁcant

digits. To circumvent such a potential cancellation, it is recommended in [10] that

the summation is carried out by coordinates, this then leads to the following tensor

product reformulation of Qq

d(f):

(2.9) Qq

d(f) =

q−1

l1=1

γq

l2=1 ···

γq

d−1

ld=1

nl1

j1=1 ···

nld

jd=1

wl1

j1···wld

jdf(xl1

j1,··· , xld

jd),

where the upper limits are deﬁned recursively as

γq

0:= q, γq

j:= γq

j−1−lj, j = 1,2,··· , d.

A FAST SG-MDI ALGORITHM FOR HIGH-DIMENSIONAL INTEGRATIPON 5

Tensor Product

3 1

1 1

2 1

Sparse Grid

3 1

Fig. 2.1: Construction of a nested sparse grid in two dimensions.

In the nested mesh case (i.e., Γl

1is nested), non-nested integration points are selected

to form the sparse grid. We remark that diﬀerent 1-d quadrature rules in (2.5) will

lead to diﬀerent sparse grid methods in (2.6). We also note that the tensor product

reformulation (2.9) will play a crucial role later in the construction of our SG-MDI

algorithm.

Figure 2.1 demonstrates the construction of a sparse grid according to the Smolyak

rule when d= 2 and q= 4. The meshes are nested, namely, Γli

1⊂Γli+1

1. The 1-d

integration-point sequence Γl1

1(l1= 1,2,3, and nl1= 3,5,9) and Γl2

1(l2= 1,2,3,

and nl2= 3,5,9) are shown at the top and left of the ﬁgure, and the tensor product

points Γ3

1⊗Γ3

1are shown in the upper right corner. From (2.6) we see that the

sparse grid rule is a combination of the low-order tensor product rule on Γi

1⊗Γj

1with

3≤i+j≤4. The point sets of these products and the resulting sparse grid Γ4

are shown in the lower half of the ﬁgure. We notice that some points in the sparse

grid Γ4

2are repeatedly used in Γi

1⊗Γj

1with 3 ≤i+j≤4. Consequently, we would

avoid the repeated points (i.e., only using the red points in Figure 2.1) and use the

reformulation (2.9) which does not involve repetition in the summation.

2.2. Examples of sparse grid methods. As mentioned earlier, diﬀerent 1-d

quadrature rules in (2.5) lead to diﬀerent sparse grid quadrature rules in (2.9). Below

we introduce four widely used sparse grid quadrature rules which will be the focus of

this paper.

Example 1: The classical trapezoidal rule. The 1-d trapezoidal rule is

文档加载中……请稍候！
如果长时间未打开，您也可以点击刷新试试。

下载文档到电脑，查找使用更方便

10 玖币 0人已下载

立即下载

摘要：

ANEFFICIENTANDFASTSPARSEGRIDALGORITHMFORHIGH-DIMENSIONALNUMERICALINTEGRATIONHUICONGZHONG∗ANDXIAOBINGFENG†Abstract.Thispaperisconcernedwithdevelopinganefficientnumericalalgorithmforfastimplementationofthesparsegridmethodforcomputingthed-dimensionalintegralofagivenfunc-tion.Thenewalgorithm,calledtheMD...

展开>> 收起<<

AN EFFICIENT AND FAST SPARSE GRID ALGORITHM FOR HIGH-DIMENSIONAL NUMERICAL INTEGRATION HUICONG ZHONGAND XIAOBING FENG.pdf

共28页,预览5页

还剩页未读，继续阅读

声明：本站为文档C2C交易模式，即用户上传的文档直接被用户下载，本站只是中间服务平台，本站所有文档下载所得的收益归上传人(含作者)所有。玖贝云文库仅提供信息存储空间，仅对用户上传内容的表现方式做保护处理，对上载内容本身不做任何修改或编辑。若文档所含内容侵犯了您的版权或隐私，请立即通知玖贝云文库，我们立即给予删除！

AN EFFICIENT AND FAST SPARSE GRID ALGORITHM FOR HIGH-DIMENSIONAL NUMERICAL INTEGRATION HUICONG ZHONGAND XIAOBING FENG

相关推荐

开通VIP享超值会员特权

作者详情

相关内容

热门标签

举报选择: