c
TU Delft, The Netherlands, 2006
MULTIGRID METHOD FOR HIGH DIMENSIONAL
ELLIPTIC EQUATIONS
H. bin Zubair∗ and C.W. Oosterlee† ∗,†Delft University of Technology, Faculty EWI
Mekelweg 4, 2628 CD, The Netherlands e-mail∗: hisham.zubair@ewi.tudelft.nl
e-mail†: oosterle@nw.twi.tudelft.nl
Key words: High-dimensions, d-multigrid, coarsening-strategies, point-smoothing Abstract. This paper discusses multigrid for high dimensional partial differential equa-tions (PDEs). We present partial grid-coarsening strategies for excellent multigrid con-vergence in the context of elliptic PDEs. We show that the multigrid concon-vergence rate can satisfactorily be brought down with the grid strategies proposed herein, coupled with weighted point smoothing schemes. A computer implementation of local Fourier smooth-ing analysis is employed to compute the optimal relaxation parameters for ω-RB Jacobi in a d-dimensional setting. The use of these relaxation parameters in the smoothing process, improves convergence in higher d. We support this by the numerical experiments in the last section, demonstrating the convergence results that we get with this proposal.
1 INTRODUCTION
Multigrid ranks among one of the best known methods for the numerical solution of partial differential equations mapped to a discrete grid [1, 2, 3, 4, 5]. Problems in ap-plied sciences, -stemming from physical systems dependent on a number of independent variables- are nowadays sometimes modelled by high dimensional partial differential equa-tions [6, 7, 8]. It has been shown [9, 10] that the Black-Scholes pricing equation for an option dependent on a basket of d assets can be reduced to the standard heat equation in d dimensions. This growth in the dimensionality of the problem renders many efficient algorithms impractical due to the asymptotic rise in the number of unknowns. A way around this so-called curse of dimensionality is the sparse grid technique. Sparse grids are always non-equidistant and so efficient solution methods for non-equidistant grids are very important in this context. Our main emphasis in this paper is on multigrid techniques through which we can overcome the deterioration of the multigrid convergence rate for values of d higher than 2 and 3.
schemes in combination with partial coarsening strategies. These strategies are based on partial doubling (h → 2h) and partial quadrupling (h → 4h) grid transfers. The proposal is to coarsen the grid continuously along the dimension(s) having the largest number of grid points, thereby reducing the stretch of the grid, till it becomes equidistant; then, to resort to full coarsening until a coarsest possible level is achieved, where an exact solution is feasible. These coarsening strategies are a mixture of partial and full transfers based on doubling and quadrupling. When these strategies are coupled with point ω-RB Jacobi, we get a nicely convergent multigrid method. Henceforth, we abbreviate multigrid for d-dimensional partial differential equations as d-multigrid.
2 THE MODEL PROBLEM AND ITS DESCRITIZATION
For analysis and experimentation we choose the d -dimensional stationary diffusion equation, with Dirichlet boundary conditions, to serve as our model problem. In what follows x is a d -tuple x = (x1, x2, · · · , xd). The continuous problem reads:
−Lu(x) = − ∂ 2 ∂x2 1 + ∂ 2 ∂x2 2 + · · · + ∂ 2 ∂x2 d u(x) = fΩ(x); x∈ Ω = (0, 1)d⊂ Rd, (1) u(x) = fΓ(x); x∈ Γ = ∂Ω; (x i ∈ {0, 1});
Subsequently, the discrete counterpart reads:
−Lhuh(x) = fhΩ(x); x∈ Ωh = (0, 1)d⊂ Z+d; (2)
uh(x) = fhΓ(x), f or x∈ Γh = ∂Ωh; (xi ∈ {0, 1})
We choose the following fourth order discretization ((4d + 1) star long-stencil) for the Laplacian. The number of cells -in the discretization grid- along the ith dimension
-represented by Ni- need not be equal to the number of cells along (say) the jthdimension.
So with hi, the mesh size along the ith dimension, the 1d stencil reads :
∂2 ∂x2 i h , 1 12h2 i [−1 16 − 30 16 − 1] + O(h4)
This discretization leads to the following matrix equation:
Ahuh = bh; (3)
3 d-MULTIGRID BASED ON POINT SMOOTHING
The core of this work is the proposed coarsening strategies which are based on a mix-ture of standard coarsening and quadrupling and which lead to efficient point smoothing based multigrid method. For anisotropic problems it is a choice to keep the point smooth-ing method but to coarsen only along a sub-set of the dimensions, precisely those that are strongly coupled. This ensures that coarsening takes place only where the errors are smooth. For isotropic problems the best strategy is to couple the point smoothing method with standard full coarsening and to use the optimal relaxation parameters for d dimensions.
Like basic multigrid for two and three dimensional problems, d -multigrid also consists of the essential components, the smoothing method and the coarse grid correction. The well known algorithm of multigrid as presented in [1] does not change for the higher dimensional case, however, the components have to be generalized to match this new situation. General multigrid algorithms are present in the literature [4, 12, 5, 2, 1, 3].
3.1 The Relaxation Method
Of the many possible choices of point smoothing based relaxation methods, we choose the ω-Red-Black Jacobi method, and in this section we give its d-dimensional adaptation. ω-RB Jacobi consists of two partial steps, each an ω-Jacobi sweep; the first one applying-to and updating only the red (odd) points, and the second one applying-applying-to and updating only the black (even) points in the grid.
This procedure indicates the need of an efficient and cheap process by which the grid G can be partitioned into the red part (Gred), and the black part (Gblk). Recall from the previous section that our grid-point-enumeration is such that the points are arranged lexicographically. We count them out in a column vector. The un-equal number of cells along different dimensions of the grid makes this partitioning process non-trivial. In this section we show how this can be accomplished by the use of the presented formulae based on Kronecker-tensor sums.
We perform this partitioning with minimal manipulation of the grid-points so that Gred
has the values at the black points ejected out and replaced with zeros and symmetrically, Gred has the values at the red positions ejected out and replaced by zeros. This ejection-process is actually where our injection operators (henceforth called ejectors) fit in.
We define the non-equidistant grid G as a vector containing the number of grid cells along each dimension, and ηiare the interior point counter along each dimension. So that:
by Er, and Eb, with: Er = d M i=1 η(d−i+1) mod 2; Eb = (Er+ 1) mod 2
L is the cummulative tensor sum of ηi, which counts the interior points along all the
dimensions , i has the reverse order (from d to 1) to match the lexicographic layout of the grid points. Due to space limitations we can only provide a 2d example, although this formulation is true in general for an abstract higher dimension d.
Example: Consider a 2d grid G = [4, 5]. In all we have 12 interior points which when counted in the lexicographic order, appear as follows:
u = [u11 u12 u13 u14 u21 u22 u23 u24 u31 u32 u33 u34]T
Evidently, G, the collection of all the points, has the following partition: Gred = {u11, u13, u22, u24, u31, u33}
Gblk = {u12, u14, u21, u23, u32, u34} Therefore according to our scheme:
ured = [ u11 0 u13 0 0 u22 0 u24 u31 0 u33 0 ]T
ublk = [ 0 u12 0 u14 u21 0 u23 0 0 u32 0 u34 ]T
and ηi in this case would be:
η1 = [1 2 3 4]T & η2 = [1 2 3]T
which leads to:
Er = (η2 ⊕ η1) mod 2
= [2 3 4 5 3 4 5 6 4 5 6 7]T mod 2
= [0 1 0 1 1 0 1 0 0 1 0 1]T
∴ Eb = [1 0 1 0 0 1 0 1 1 0 1 0]T
(5)
These red and black point ejectors now can be used to partition the grid as described. Once this partition is obtained, ω-RB Jacobi relaxation sweeps are trivial to perform, in that, they are no different than their 2d counterparts.
We also employ optimal relaxation parameters ωopt in the relaxation process for
3.2 Coarsening strategies to handle anisotropies
We present two grid-adaptive coarsening procedures as our test cases, both of which have shown excellent convergence results. We call them strategy-1 and strategy-2. In Strategy-1 we first identify the dimensions having the largest number of grid points and then we coarsen simultaneously along these dimensions through an (h → 2h) transfer. When this strategy is employed recursively at each level, it makes the grid equidistant after a few simultaneous partial coarsenings. At this stage we start full coarsening (h → 2h transfers along all dimensions) until we finally reach the coarsest possible level and then we solve exactly. In the last section we present experimental results for multigrid convergence that evaluate these strategies.
Strategy-2 is similar to Strategy-1 in that it also employs simultaneous partial coarsen-ing, however, here we apply quadrupling in the partial phase and standard coarsening in the full phase (after achieving an equidistant grid). This strategy gives good convergence when employed in conjunction with optimal relaxation parameters, and, is cheaper than strategy-1 because of quadrupling in the partial phase. We also point out that a strategy based on full quadrupling should not be employed in general dimensions (d > 2) as it results in poor convergence, even with the use of optimal relaxation parameters.
3.3 Coarse-grid discretization
An important component in the coarse grid correction process is the choice of the coarse-grid operator LH. In this paper we use the coarse-grid analog of the discrete
operator on the fine-grid. Once the next coarser-grid is decided we discretize the Laplacian using the same discrete stencils as presented in section 2.
3.4 The transfer Operators
We employ the d-dimensional adaptations of the Full-Weighting (FW) restriction op-erator, and of the bilinear interpolation operator in 2d; for the intergrid transfers of the grid functions. In this section we present a tensor formulation to generate the restric-tion and prolongarestric-tion operator matrices. For completeness we first menrestric-tion [1] that a 2d FW restriction operator is the Kronecker tensor product of the x1 and x2 directional
1-dimensional FW operators. Likewise for a general FW restriction operator matrix, we have a formula based on Kronecker tensor products for building up a FW restriction operator matrix R, which reads:.
We now define the quantities involved in (6) for the dummy subscript a. Ia is the identity matrix of order (a − 1) × (a − 1).
Oa is the 1d FW restriction operator matrix, order = a2 − 1 × (a − 1).
T= [k1, k2, · · · , kd] is the coarsening request, ki is the count of (h → 2h) transfers along
the ith dimension. For reasons of space it is not possible to verify (in this manuscript)
eqn. (6) for any realistic example, however it is trivial to verify the same with a matrix manipulation software package.
Once the FW restriction operator matrix in d-dimensions is set, the prolongation (d-linear interpolation) operator matrix can be obtained by the following relation:
P= 2d(RT) (7)
Note that the FW restriction operator given by (6) provides the required operator ma-trix for any number of coarsenings along any number of dimensions for an abstract d-dimensional problem.
3.5 Computational Work for d-multigrid
The practical feasibility of a d-multigrid method also has to take into account an es-timate of the computational work that it involves. It is important to point out that the coarsening strategies that we use do not employ the same transfers at all grid levels. The strategies are hybrid in the sense that they are composed of a mixture of standard and quadrupling transfers, moreover, they are grid-adaptive, implying that the coarsening de-cision is made depending on the existing stretch of the grid. Note that in our coarsening strategies the number of dimensions along which the grid is coarsened, increases (or re-mains constant in the case of an equidistant grid) at each successive level. This gives the following [11], Wl 6 2js (2js− γ)CMl Wl 6 τ (τ − γ)CMl; (τ = 2 js) (8)
Here j represents the number of dimensions coarsened and s = 1 for (h → 2h) transfers and s = 2 for quadrupling. The worst case would be when the grid is highly stretched along a single dimension, which implies that coarsening takes place only along j = 1 dimension. With quadrupling, i.e. s = 2 this still renders τ > 2, which implies that the complexity of the method is still O(M) even with W cycles. This feature makes quadrupling particularly attractive in higher dimensions for unequidistant grids.
If the problem is isotropic i.e. the grid is equidistant, we essentially have s = 1 and j = d, which gives the work estimates in Table 1,
γ d = 2 d = 3 d = 4 d = 5 d = 6 1 43 CMl 87 CMl 1516 CMl 3231 CMl 6463 CMl
2 2 CMl 43 CMl 87 CMl 1615 CMl 3231 CMl
4 O(Mllog2Ml) 2 CMl 34 CMl 87 CMl 1615 CMl
Table 1: Work estimates for d-multigrid
worst case (grid highly stretched along a single dimension), Strategy 2 still results in an O(M) multigrid method.
4 OPTIMAL RELAXATION PARAMETERS FOR ω-RB JACOBI
A local Fourier smoothing analysis of ω-RB Jacobi for 2 and 3 dimensional problems is present in the literature [4, 3]. In [11] we have shown how to extend it to a general d dimensional setting. From a computer implementation of this analysis we obtain relax-ation parameters which for an assortment of cases are provided in Table 3. They enhance convergence considerably in dimensions higher than 2.
Doubling (h → 2h) d µ(ω) = µ(1) ωopt µ(ωopt) 2 0.284 1.064 0.169 3 0.460 1.135 0.210 4 0.571 1.185 0.238 5 0.646 1.224 0.262 6 0.698 1.256 0.283 Quadrupling (h → 4h) µ(ω) = µ(1) ωopt µ(ωopt) 0.761 1.312 0.399 0.835 1.378 0.470 0.874 1.424 0.519 0.898 1.458 0.555 0.915 1.485 0.584
Table 2: ωopt and µ, equidistant grid, O(h4) long-stencil, d represents the number of dimensions
Doubling (h → 2h) d µ(ω) = µ(1) ωopt µ(ωopt) 2 0.174 1.009 0.166 3 0.267 1.056 0.219 4 0.345 1.120 0.249 5 0.409 1.180 0.276 6 0.462 1.230 0.304 Quadrupling (h → 4h) µ(ω) = µ(1) ωopt µ(ωopt) 0.640 1.234 0.315 0.692 1.255 0.387 0.730 1.283 0.411 0.761 1.311 0.402 0.785 1.333 0.397
Table 3: ωoptand µ, 128 × 32(d−1) grid, O(h4) long-stencil, d represents the number of dimensions
Experiment 2: d = 2; G = 1282 V W V W [µ(1)]2 = 0.08 [µ(1.06)]2 = 0.03 O(h4) 0.13/9 0.09/8 0.12/8 0.07/7 Experiment 2: d = 3; G2 = 1283, G4 = 643 V W V W [µ(1)]2 = 0.21 [µ(1.13)]2 = 0.21 O(h4) 0.26/12 0.22/11 0.14/9 0.08/8 Experiment 3: d = 4; G2 = 644, G4 = 324 V W V W [µ(1)]2 = 0.33 [µ(0.24)]2 = 0.06 O(h4) 0.39/16 0.34/14 0.18/10 0.09/8 Experiment 4: d = 5; G = 165 V W V W [µ(1)]2 = 0.41 [µ(1.24)]2 = 0.09 O(h2) ρ(M h)/n 0.38/16 0.38/15 0.19/10 0.09/8 Experiment 5: d = 6; G = 86 V W V W [µ(1)]2 = 0.48 [µ(1.28)]2 = 0.12 O(h2) ρ(M h)/n 0.35/15 0.34/15 0.10/9 0.10/8
Table 4: Result-summary of numerical experiments on equidistant grids.
5 NUMERICAL EXPERIMENTS
Now we present the results of some of our numerical experiments. These results are drawn both for Strategy-1 and Strategy-2. The multigrid convergence factor ρ(Mh) [1],
displayed in Tables 4-9 is presented against the number of cycles that it took to converge to a tolerance value of 10−6. Results are included for both the V and the W cycles. G
represents the discretization grid, each component stands for the number of cells along that dimension. Results are included for the fourth order long stencil upto the fourth dimension and then for the second order stencil for d = 5 & d = 6.
Table 4 presents experimental results for equidistant grids. Tables 5-9 present results for different non-equidistant grids, some of them stretched along a single dimension, thus simulating a worst case in the context of computational expense. The multigrid method remains O(M) with strategy 2 (partial quadrupling) even in the worst cases, furthermore the results show that optimal parameters positively influence the multigrid convergence factor as well as the total number of multigrid cycles, with the anisotropic case as well.
6 CONCLUSIONS
Strategy 1 G = [512, 32] V W V W [µ(1)]2 = 0.01 [µ(0.99)]2 = 0.01 O(h4) 0.10/8 0.03/6 0.11/8 0.04/6 Strategy 2 G = [512, 32] V W V W [µ(1)]2 = 0.35 [µ(1.22)]2 = 0.05 0.34/14 0.31/13 0.15/9 0.09/8
Table 5: Observed ρ(Mh) against the number of iterations, for 2d
Strategy 1 G = [128, 32, 32] V W V W [µ(1)]2 = 0.07 [µ(1.06)]2 = 0.05 O(h4) 0.24/11 0.04/6 0.16/9 0.06/7 Strategy 2 G = [128, 32, 32] V W V W [µ(1)]2 = 0.48 [µ(1.25)]2 = 0.15 0.34/14 0.35/14 0.17/9 0.16/8
Table 6: Observed ρ(Mh) against the number of iterations, for 3d
Strategy 1 G = [128, 32, 32, 32] V W V W [µ(1)]2 = 0.10 [µ(1.08)]2 = 0.03 O(h4) 0.33/13 0.04/6 0.20/10 0.07/7 Strategy 2 G = [128, 32, 32, 32] V W V W [µ(1)]2 = 0.53 [µ(1.28)]2 = 0.17 0.38/16 0.37/15 0.20/10 0.16/8
Table 7: Observed ρ(Mh) against the number of iterations, for 4d
Strategy 1 G = [128, 8, 8, 32, 8] V W V W [µ(1)]2 = 0.04 [µ(1.04)]2 = 0.02 O(h2) 0.28/13 0.03/5 0.14/9 0.04/6 Strategy 2 G = [128, 8, 8, 32, 8] V W V W [µ(1)]2 = 0.40 [µ(1.24)]2 = 0.09 0.54/22 0.27/12 0.26/12 0.10/7
Table 8: Observed ρ(Mh) against the number of iterations, for 5d
Strategy 1 Experiment 10a: d = 6; G = 1282× 44 V W V W [µ(1)]2 = 0.08 [µ(1.06)]2 = 0.03 O(h2) 1.11/9 0.06/6 0.10/8 0.05/6 Strategy 2 Experiment 10b: d = 6; G = 1282× 44 V W V W [µ(1)]2 = 0.52 [µ(1.33)]2 = 0.11 0.52/22 0.52/20 0.26/11 0.11/8
relaxation parameters, and suitable coarsening steps we render nicely convergent multi-grid methods. Point ω-RB Jacobi proves to be an excellent relaxation method. The relaxation parameters considerably enhance the multigrid convergence for d > 2. Partial quadrupling gives an O(M) W -cycle multigrid method even in the worst case of partial coarsening solely in one direction.
REFERENCES
[1] Trottenberg U., Oosterlee C.W., and Sch¨uller A. Multigrid. Academic Press, 2001. [2] Wesseling P. An Introduction to Multigrid Methods. Pure and Applied Mathematics.
John Wiley & Sons, 1992.
[3] Wienands R. and Wolfgang J. Practical Fourier Analysis for Multigrid Methods. Numerical Insights. Chapman & Hall/CRC, 2005.
[4] St¨uben K. and Trottenberg U. On the construction of fast solvers for elliptic equa-tions. Technical report, von Karman Institute for Fluid Dynamics, Rhode-Saint-Genese, Belgium, 1982.
[5] Thole C.A. and Trottenberg U. Basic smoothing procedures for the multigrid treat-ment of elliptic 3-d operators. Appl. Math. Comp., 19:pp 333–345, 1986.
[6] Elf J., L¨otstedt P., and Sj¨oberg P. Problems of high dimension in molecular biology. In Proceedings of the 19th GAMM-Seminar Leipzig, pages 21–30, 2003.
[7] Yserentant H. Sparse grid spaces for the numerical solution of the electronic schr¨odinger equation. Numer. Math., pages pp 381–389, 2005.
[8] Reisinger C. and Wittum G. On multigrid for anisotropic equations and variational inequalities. Computing and Visualization in Science, 2004.
[9] Kwok Y. K. Mathematical models of financial derivatives. Springer Finance. Springer-Verlag Singapore, Singapore, second edition, 1998.
[10] You1an Z., Xiaonan W., and I-Liang C. Derivative Securities and Difference Methods. Springer Finance. Springer Science+Business Media Inc. USA, 2004.
[11] bin Zubair H., Wienands R., and Oosterlee C.W. Multigrid for high dimensional elliptic partial differential equations on non-equidistant grids. Technical report, Delft University of Technology, 2006.
