# Local Fourier analysis of Balancing Domain Decomposition by Constraints algorithms

Local Fourier analysis is a commonly used tool for the analysis of multigrid and other multilevel algorithms, providing both insight into observed convergence rates and predictive analysis of the performance of many algorithms. In this paper, for the first time, we adapt local Fourier analysis to examine variants of two- and three-level balancing domain decomposition by constraints (BDDC) algorithms, to better understand the eigenvalue distributions and condition number bounds on these preconditioned operators. This adaptation is based on a modified choice of basis for the space of Fourier harmonics that greatly simplifies the application of local Fourier analysis in this setting. The local Fourier analysis is validated by considering the two dimensional Laplacian and predicting the condition numbers of the preconditioned operators with different sizes of subdomains. Several variants are analyzed, showing the two- and three-level performance of the "lumped" variant can be greatly improved when used in multiplicative combination with a weighted diagonal scaling preconditioner, with weight optimized through the use of LFA.

## Authors

• 12 publications
• 22 publications
• 10 publications
01/03/2020

### Tuning Multigrid Methods with Robust Optimization

Local Fourier analysis is a useful tool for predicting and analyzing the...
12/19/2019

### On the stability of Scott-Zhang type operators and application to multilevel preconditioning in fractional diffusion

We provide an endpoint stability result for Scott-Zhang type operators i...
01/21/2020

### Balancing domain decomposition by constraints associated with subobjects

A simple variant of the BDDC preconditioner in which constraints are imp...
08/26/2019

### A local Fourier analysis of additive Vanka relaxation for the Stokes equations

Multigrid methods are popular solution algorithms for many discretized P...
06/10/2022

### Linear peridynamics Fourier multipliers and eigenvalues

A characterization for the Fourier multipliers and eigenvalues of linear...
01/14/2019

### A Multilevel Approach for the Performance Analysis of Parallel Algorithms

We provide a multilevel approach for analysing performances of parallel ...
12/06/2020

### Fourier-domain Variational Formulation and Its Well-posedness for Supervised Learning

A supervised learning problem is to find a function in a hypothesis func...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

Domain decomposition methods are well-studied approaches for the numerical solution of partial differential equations both experimentally and theoretically

[5, 13, 15, 36], due to their efficiency and robustness for many large-scale problems, and the need for parallel algorithms. Among the many families of domain decomposition algorithms are Neumann-Neumann [36], FETI-DP [16], Schwarz [15, 36], and Optimized Schwarz [13, 20]. Balancing domain decomposition by constraints (BDDC) is one family of non-overlapping domain decomposition methods. While BDDC was first introduced by Dohrmann in [10], several variants have recently been proposed. BDDC-like methods have been successfully applied to many PDEs, including elliptic problems [25], the incompressible Stokes equations [24, 26], H(curl) problems [12], flow in porous media [38], and the incompressible elasticity problem [11, 32]. Theoretical analysis of BDDC has primarily been based on finite-element approximation theory [8, 14, 29, 30]. It has been shown that the condition number of the preconditioned BDDC operator can be bounded by a function of (where is the meshsize, and is the subdomain size), independent of the number of subdomains [29, 30]. A nonoverlapping domain decomposition method for discontinuous Galerkin based on the BDDC algorithm is presented in [7]

, and the condition number of the preconditioned system is shown to be bounded by similar estimates as those for conforming finite element methods. BDDC methods in three- or multilevel forms have also been developed

[31, 39, 40], and good implementations are available, for example, [1, 34, 43].

Since BDDC algorithms are widely used to solve many problems with high efficiency and parallelism, better understanding of how this methodology works is useful in the design of new algorithms. Local Fourier analysis (LFA), first introduced by Brandt [6] and well-studied for multigrid methods [9, 35, 37, 41, 42], is an analysis framework that provides predictive performance estimates for many multilevel iterations and preconditioners. Early application of LFA was mainly focused on scalar problems or systems of PDEs with collocated discretizations. Standard LFA [37] cannot be directly applied to higher-order and staggered-grid discretizations, since the analysis depends on Toeplitz operator structure inherited from the mesh and discretization. There is a long history of generalization of “standard” LFA to account for more general structure of the discrete operator and/or multigrid algorithm, beginning with analysis of red-black (and, later, multicolour) relaxation [22, 33, 37]. This work has been generalized in recent years, leading to the idea of “periodic stencils” [3] and similar approaches for PDEs with random coefficients [21]. From a different perspective, similar tools have arisen to account for the structure of coupled systems of PDEs, beginning with [4]. This work was expanded by MacLachlan and Oosterlee [28], to account for the structure of overlapping relaxation schemes for the Stokes Equations. In a third setting, the mode analysis of certain space-time multigrid methods, Friedhoff and MacLachlan again use a similar approach to handle coarsening in time [18, 19]. In this work, we show how a generalization of these approaches can be applied to the structure of domain decomposition algorithms.

To our knowledge, there has been no research applying local Fourier analysis to BDDC-like algorithms. The same is true of the closely related finite element tearing and interconnecting dual-primal (FETI-DP) methodology [16, 17, 23]. Here, we adapt LFA to this domain decomposition method by borrowing tools from LFA for systems of PDEs and other block-structured settings. Noting that the stencil for the domain decomposition method depends on the size of the subdomains, an adaptation of the standard basis is useful, which we present in Section 4.1

. Because LFA can reflect both the distribution of eigenvalues and associated eigenvectors of a preconditioned operator, here, we adopt LFA to analyze variants of the common “lumped” and “Dirichlet” BDDC algorithms, based on

[27], to guide construction of these methods. To do this, we use a modified basis as in [3, 21, 28] for the Fourier analysis that is well-suited for application to domain decomposition preconditioners.

Applying the two-level BDDC algorithm requires the solution of a Schur complement equation (coarse problem), which usually poses some difficulty with increasing problem size. Two- and three-level variants are, thus, considered in this paper. However, as is well-known in the literature, bounds on the performance of BDDC degrade sharply from two-level to three-level methods, particularly for large values of . Since our analysis shows that the largest eigenvalues of the preconditioned operator for the lumped BDDC algorithm are associated with oscillatory modes, we propose variants of standard BDDC based on multiplicative preconditioning and multigrid ideas. From the condition numbers offered by LFA, we can easily compare the efficiency of these variants. Furthermore, LFA can provide optimal parameters for these multiplicative methods, helping tune and understand sensitivity to the parameter choice.

This paper is organized as follows. In Section 2, we introduce the finite element discretization of the Laplace problem in two dimensions and the lumped and Dirichlet preconditioners. Two- and three-level preconditioned operators are developed in Section 3. In Section 4, we discuss the Fourier representation of the preconditioned operators. Section 5 reports LFA-predicted condition numbers of the BDDC variants considered here. Conclusions are presented in Section 6.

## 2 Discretization

We consider the two-dimensional Laplace problem in weak form: Find such that

 a(u,v)=∫Ω∇u⋅∇vdΩ=⟨f,v⟩,∀v∈V, (1)

where is a bounded domain with Lipschitz boundary . Here, we consider the Ritz-Galerkin approximation over , the space of piecewise bilinear functions on a uniform rectangular mesh of . The corresponding linear system of equations is given as

 Ax=b. (2)

We partition the domain, , into nonoverlapping subdomains, , where each subdomain is a union of shape regular elements and the nodes on the boundaries of neighboring subdomains match across the interface . The interface of subdomain is defined by . Here, we consider with both a discretization mesh (with meshsize ) and subdomain mesh (with meshsize ) given by uniform grids with square elements or subdomains.

The finite-element space can be rewritten as where is the product of the subdomain interior variable spaces . Functions in are supported in the subdomain and vanish on the subdomain interface . is the space of traces on of functions in . Then, we can write the subdomain problem with Neumann boundary conditions on as

 A(i)x(i)=⎛⎝A(i)IIA(i)TΓIA(i)ΓIA(i)ΓΓ⎞⎠(x(i)Ix(i)Γ)=(b(i)Ib(i)Γ), (3)

where , and is the (conjugate) transpose. Then, the global problem (2) can be assembled from the subdomain problems (3) as

 A=N∑i=1R(i)TA(i)R(i),andb=N∑i=1R(i)Tb(i),

where

is the restriction operator from a global vector to a subdomain vector on

.

### 2.1 A Partially Subassembled Problem

In order to describe variants of the BDDC methods, we first introduce a partially subassembled problem, following [27], and the corresponding space of partially subassembled variables,

 ^Vh=VΠ,h⨁Vr,h, (4)

where

is spanned by the subdomain vertex nodal basis functions (the coarse degrees of freedom). The complementary space,

, is the product of the subdomain spaces , which correspond to the subdomain interior and interface degrees of freedom and are spanned by the basis functions which vanish at the coarse-grid degrees of freedom. For a mesh, the degrees of freedom in are those corresponding to the circled nodes at the left of Figure 1, while the degrees of freedom in correspond to all interior nodes, plus duplicated (broken) degrees of freedom along subdomain boundaries.

The partially subassembled problem matrix, corresponding to the variables in the space , is obtained by assembling the subdomain matrices (3) only with respect to the coarse-level variables; that is,

 ^A=N∑i=1(¯R(i))TA(i)¯R(i), (5)

where is a restriction from space to . The central idea behind BDDC preconditioners is to use as an approximation to that is more readily inverted, since the subdomain problems can be “condensed” to a solve on the underlying coarse grid, as described below.

###### Remark 1

The above “broken” space can, in fact, be viewed by considering the discretization of the PDE independently on each subdomain and only “assembling” the DOFs at the corners of the subdomains. Then, looking at this “subassembled” problem, the interface between any two subdomains has independent DOFs associated with each subdomain at the same locations, which are normally assembled to formulate the global system. This is different than the standard overlapping domain decomposition, which uses only one DOF at these locations in the overlap between subdomains and contributes this value to each of the shared subdomains. While can be readily defined in terms of , the reverse is not true.

### 2.2 Lumped and Dirichlet Preconditioners

In order to define the preconditioners under consideration for (2), we introduce a positive scaling factor, , for each node on the interface of subdomain Let be the set of indices of the subdomains that have on their boundaries. Define , where is the cardinality of . The scaled injection operator, , is defined so that each column of corresponds to a degree of freedom of the global problem (2). For subdomain interior and coarse-level variables, the corresponding column of has a single entry with value 1. Columns that correspond to an interface degree of freedom (the set of nodes in ) have non-zero entries each of .

Based on the partially subassembled problem, the first preconditioner introduced for solving (2) is

 M−11=RT1^A−1R1.

The preconditioned operator has the same eigenvalues as the preconditioned FETI-DP operator with a lumped preconditioner, except for some eigenvalues equal to 0 and 1 [17, 27]. We refer to as the lumped preconditioner. Note that corresponds to a naive approximation of by , using equal weighting of residual and corrections to/from the duplicated degrees of freedom in the partially broken space.

A similar preconditioner for

augments this by using discrete harmonic extensions in the restriction and interpolation operators

[27], giving

 M−12=(RT1−HJD)^A−1(R1−JTDHT):=R2, (6)

where is the direct sum of local operators , which map the jump over a subdomain interface (given by ) to the interior of the subdomain by solving a local Dirichlet problem, and gives zero for other values. For any given , the component of on subdomain is given by

 (JTDv(x))(i)=∑j∈Nx(δj(x)v(i)(x)−δi(x)v(j)(x)),∀x∈Γi. (7)

Extending the interface values using the discrete harmonic extension minimizes the energy norm of the resulting vector [36], giving a better stability bound. Furthermore, the preconditioned operator has the same eigenvalues as the BDDC operator [25], except for some eigenvalues equal to 1 [27]. We refer to as the Dirichlet preconditioner. Note that differs from only by local operators in the operators used to map between the fully assembled and partially broken spaces.

Standard bounds (see, e.g., [27]) on the condition numbers of the preconditioned operators are that, for , there exists such that and, for , there exists such that .

###### Remark 2

While the preconditioners given by have some notational similarity to Galerkin corrections, it is important to note that acts on a higher-dimensional space than , and does not satisfy for either “restriction” operator. The essential difference between the two preconditioners is in how “glues” the solution over the broken space back into the usual continuous space. does this through a simple partition of unity, while minimizes subdomain energy when propagating mismatches along the subdomain boundaries.

## 3 Two- and Three-level Variants

In both of the above preconditioned operators, we need to solve the partially subassembled problems, that we write in block form

 (Arr^ATΠr^AΠrAΠΠ)^A(^xr^xΠ)^x=(Arr0^AΠr^SΠ)(IA−1rr^ATΠr0I)(^xr^xΠ)=(^dr^dΠ)^d, (8)

where

is the identity matrix,

and contains the subdomain interior and interface degrees of freedom, and corresponds to the coarse-level degrees of freedom, which are located at the corners of the subdomains. We write in (8) in factorization form to easily separate the action on the coarse degrees of freedom, and to find the corresponding symbol of . If we define

 P=(−A−1rr^ATΠrI),

then the Schur complement is simply the Galerkin coarse operator, , and the block-factorization solve for is equivalent to a two-level additive multigrid method with exact -relaxation using

 SF=(A−1rr000),

to define an “ideal” relaxation scheme to complement the coarse-grid correction defined by .

In the partially subassembled problem (8), we need to solve a coarse problem related to . We can either solve this coarse problem exactly (corresponding to a two-level method, where the Schur complement is inverted exactly) or inexactly (as a three-level method), where the lumped and Dirichlet preconditioners defined above are used recursively to solve this problem.

### 3.1 Exact and Inexact Solve for the Schur Complement

Let

 KLD=(Arr0^AΠr^SΠ), KU=(IA−1rr^ATΠr0I),

and note that and are formed from the LDU factorization of as in (8). For , let denote the preconditioned operators for two- and three-level variants of BDDC, where and denote using and (with as preconditioners for the fine and coarse problems, respectively, where stands for applying the preconditioner to the Schur complement matrix, . By standard calculation, we can write

 Gi,j=RTiK−1UPjK−1LDRiA,

with

 Pj=(I00M−1s,j^SΠ).
###### Remark 3

When , is a two-level method, solving the Schur complement problem exactly, as . Note that, for the three-level variants ,

 PjK−1LD=(I00M−1s,j^SΠ)(A−1rr0−^S−1Π^AΠrA−1rr^S−1Π)=(A−1rr0−M−1s,j^AΠrA−1rrM−1s,j),

corresponding to replacing in (8) by representing the inexact solve. From this, it is clear that can be applied without directly applying the inverse of .

###### Theorem 1

The eigenvalues of are real and bounded below by 1.

###### Proof

This result is well-known [8, 25] for the two-level methods, . For , note that

 KLDP−1jKU = (Arr0^AΠ^SΠ)(I00^S−1ΠMs,j)(IA−1rrATΠ,r0I) = KTU(Arr00Ms,j)KU,

which is symmetric and positive definite (SPD), ensuring the eigenvalues of are real. Similarly, can be rewritten as

 ^A = (I0^AΠrA−1rrI)(Arr00^SΠ)(IA−1rr^ATΠr0I)=KTU(Arr00^SΠ)KU.

Since is SPD, so is . We form from in the same way. As for the two-level case [8, 25], we can bound

 uTrMs,jur≤uTr^SΠur,∀ur.

Writing , we have

 uTIArruI+uTrMs,jur≤uTIArruI+uTr^SΠur,

 uT(Arr00Ms,j)u≤uT(Arr00^SΠ)u. (9)

Now, for any , taking in (9), we have

which is

 yTKLDP−1jKUy≤yT^Ay.

According to [27, Theorem 4], we then have

 ⟨u,u⟩A≤⟨Gi,ju,u⟩A,∀u,

which means that the smallest eigenvalue of is not less than 1.

###### Remark 4

In [27], represents one or several multigrid V-cycles or W-cycles for solving the global partially subassembled problem (8). In our setting, we use the lumped and Dirichlet preconditioners recursively, which also keeps the lower bound, 1, for the inexact solve.

###### Remark 5

Standard bounds (see, e.g., [27, 40]) on the condition numbers of the three-level preconditioned operators are that there exists such that , where and .

### 3.2 Multiplicative Preconditioners

As we shall see, the bounds above are relatively sharp and the performance of both preconditioners degrades with subdomain size and number of levels. To attempt to counteract this, we consider multiplicative combinations of these preconditioners with a simple diagonal scaling operator, mimicking the use of weighted Jacobi relaxation in classical multigrid methods. We use to denote the multiplicative preconditioned operator based on with diagonal scaling on the fine level. Here,

 Gfi,j=Gi,j+ωD−1A(I−Gi,j),i=1,2,j=0,1,2, (10)

where is the diagonal of and is a chosen relaxation parameter. Note that , so represents the multiplicative combination given.

###### Theorem 2

The eigenvalues of are real.

###### Proof

First, recall that is SPD, so it has a unique SPD matrix square root, . Note that . From Theorem 1, we know that is symmetric, so is symmetric as well.

is similar to , which can be rewritten as

 A12(Gfi,j−I)A−12 = A12(I−ωD−1A)(Gi,j−I)A−12 = (I−ωA12D−1A12)(A12Gi,jA−12−I)

From Theorem 1, we know that the eigenvalues of are not less than 1. Thus, is symmetric positive semi-definite, and has a unique symmetric positive semi-definite square root. Using the fact that , we have

 λ(Gfi,j−I) = λ((I−ωA12D−1A12)(A12Gi,jA−12−I)) = λ((A12Gi,jA−12−I)12(I−ωA12D−1A12)(A12Gi,jA−12−I)12)

Note that is symmetric. Thus, the eigenvalues of are real and so are those of .

Another variant is the use of multiplicative preconditioning on the coarse level with a similar diagonal scaling. We use to denote the resulting multiplicative preconditioner. Here,

 (11)

where

 Pcj=(I00Gc,j),

in which

 Gc,j=M−1s,j^SΠ+ωD−1s^SΠ(I−M−1s,j^SΠ),

where is the diagonal of .

Instead of using a single sweep of Jacobi in , we can consider a symmetrized Jacobi operator , where ; that is,

 Gsc,j=Gc,j+ω2(I−Gc,j)D−1s^SΠ,

then changes to

 Gs,ci,j=RTiK−1UPs,cjK−1LDRiA,i,j=1,2, (12)

for defined as but with in the -block.

Finally, we can also apply the multiplicative operators based on diagonal scaling on both the fine and coarse levels. We denote this as

 Gf,ci,j=Gci,j+ω2D−1A(I−Gci,j),i=1,2,j=1,2, (13)

where is the diagonal of and is a chosen relaxation parameter.

###### Remark 6

The eigenvalues of operators , and are not generally all real. When , is similar to a symmetric matrix and we observe a real spectrum in numerical results, but not for all cases when . In most situations, and have complex eigenvalues.

In the following, we focus on analyzing the spectral properties of the above preconditioned operators by local Fourier analysis [37]. The main focus of this work is on the operators and , because the Fourier representations of other operators are just combinations of these three and some simple additional terms.

## 4 Local Fourier Analysis

To apply LFA to the BDDC-like methods proposed here, we first review some terminology of classical LFA. We consider a two-dimensional infinite uniform grid, , with

 Gh={xi,j:=(xi,xj)=(ih,jh),(i,j)∈Z2}, (14)

and Fourier functions on , where and . Let be a Toeplitz operator acting on as

 Lh∧=[sκ]h(κ=(κ1,κ2)∈Z2);Lhwh(x)=∑κ∈Vsκwh(x+κh),

with constant coefficients , where is a function on . Here, is taken to be a finite index set. Note that since is Toeplitz, it is diagonalized by the Fourier modes .

###### Definition 1

If for all grid functions, ,

 Lhψ(θ,x)=˜Lh(θ)ψ(θ,x),

we call the symbol of .

###### Remark 7

In Definition 1, the operator acts on a single function on , so is a scalar. For an operator mapping vectors on to vectors on , the symbol will be extended to be a matrix.

### 4.1 Change of Fourier Basis

Here, we discuss LFA for a domain decomposition method. While the classical basis set for LFA, denoted below, could be used, we find it is substantially more convenient to make use of a transformed “sparse” basis, introduced here as . This basis allows a natural expression of the periodic structures in domain decomposition preconditioners that vary with subdomain size. We treat each subdomain problem as one macroelement patch, and each subdomain block in the global problem is diagonalized by a coupled set of Fourier modes introduced in the following. Because each subdomain has the same size, , we consider the high and low frequencies for coarsening by factor , given by

 θ∈Tlow=[−πp,πp)2,θ∈Thigh=[−πp,(2p−1)πp)2\[−πp,πp)2.

Let , where and for . For any given , we define the -dimensional space

 Eh(θ(0,0)):=span{ψ(θ(q,r),xs,t)=eιθ(q,r)⋅xs,t/h:q,r=0,1,⋯,p−1}, (15)

as the classical space of Fourier harmonics for factor coarsening.

For any , we consider a grid function defined as a linear combination of the basis functions for with frequencies and coefficients as

 es,t:=p−1∑q,r=0βq,rψ(θ(q,r),xs,t).

We note that any index has a unique representation as where and . From (15), we have

 epm+k,pn+ℓ = p−1∑q,r=0βq,reι(θ(0)1+2πqp)xs/heι(θ(0)2+2πrp)xt/h = p−1∑q,r=0βq,reιθ(0)1xs/heι2πq(pm+k)peιθ(0)2xt/he2πr(pn+ℓ)p = p−1∑q,r=0βq,reι2πqkpeιθ(0)1xs/he2πrℓpeιθ(0)2xt/h = (p−1∑q,r=0βq,reι2πqkpe2πrℓp)(eιθ(0,0)⋅xs,t/h).

Thus, we can write

 epm+k,pn+ℓ=^βk,ℓeιθ⋅xs,t/H, (16)

with

 θ=pθ(0,0),and^βk,ℓ=p−1∑q,r=0βq,reι2πqkpe2πrℓp. (17)

In other words, for any point with and , can be reconstructed from a single Fourier mode with coefficient . Thus, on the mesh defined in (14), the periodicity of the basis functions in can also be represented by a pointwise basis on each -block.

Based on (16), we consider a “sparse” -dimensional space as follows

 FH(θ):=span{φk,ℓ(θ,xs,t)=eιθ⋅xs,t/Hχk,ℓ(xs,t):k,ℓ=0,1,⋯,p−1}, (18)

where and

 χk,ℓ(xs,t)={1,if mod(s,p)=k,andmod(t,p)=ℓ,0,otherwise.

Note that, with this notation, (16) can be rewritten as

 epm+k,pn+ℓ=^βk,ℓφk,ℓ(θ,xs,t). (19)
###### Theorem 3

and are equivalent.

###### Proof

While the derivation above shows directly that , we revisit this calculation now to show that the mapping is invertible and, hence, as well.

Let be an arbitrary vector with size , denoted as

Then, we define a vector, , based on (17), as follows

 ^X=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝^X0^X1⋮^Xp−2^Xp−1⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠,where^Xℓ=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝^β0,ℓ^β1,ℓ⋮^βp−2,ℓ^βp−1,ℓ⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠,ℓ=0,1,⋯,p−1,

in which

 ^βk,ℓ=p−1∑r=0(p−1∑q=0βq,reι2πqkp)e2πrℓp,q,r=0,1,⋯,p−1.

Let be the matrix of this transformation, , and

Note that defines a vector whose -th entry is and, thus, we see that .

Note that is a Vandermonde matrix based on values , where . It is obvious that if . Consequently, . Thus, is invertible, and so is . It follows that and are equivalent.

###### Remark 8

Let , be the primitive -th root of unity, and note that . Thus,

is the unitary discrete Fourier transform (DFT) matrix with