DeepAI

# On condition numbers of symmetric and nonsymmetric domain decomposition methods

Using oblique projections and angles between subspaces we write condition number estimates for abstract nonsymmetric domain decomposition methods. In particular, we design and estimate the condition number of restricted additive Schwarz methods. We also obtain non-negativity of the pre-conditioner operator. Condition number estimates are not enough for the convergence of iterative method such as GMRES but these bounds may lead to further understanding of restricted methods.

10/08/2021

### Three decompositions of symmetric tensors have similar condition numbers

We relate the condition numbers of computing three decompositions of sym...
03/20/2020

### Analysis of the SORAS domain decomposition preconditioner for non-self-adjoint or indefinite problems

We analyze the convergence of the one-level overlapping domain decomposi...
04/01/2021

### An abstract theory of domain decomposition methods with coarse spaces of the GenEO family

Two-level domain decomposition methods are preconditioned Krylov solvers...
08/19/2020

### On the condition number of the total least squares problem with linear equality constraint

This paper is devoted to the condition number of the total least squares...
11/24/2021

### Sensitivity Analysis of Domain Decomposition method for 4D Variational Data Assimilation (DD-4DVAR DA)

We prove consistence, convergence and stability of the Domain Decomposit...
07/26/2019

### Learning high-dimensional additive models on the torus

In this paper we study the multivariate ANOVA decomposition for 1-period...
01/21/2020

### Balancing domain decomposition by constraints associated with subobjects

A simple variant of the BDDC preconditioner in which constraints are imp...

## 1 Introduction

The restricted additive Schwarz (RAS) was originally introduced by Cai and Sarkis in [4] in 1999. RAS outperforms the classical additive Schwarz (AS) preconditioner in the sense that it requires fewer iterations, as well as lower communication and CPU time costs when implemented on distributed memory computers, [4]. Unfortunately, RAS in its original form is nonsymmetric, and therefore the conjugate gradient (CG) method cannot be used. Pursuing the analysis of RAS, several interesting methods have been developed. Some of these versions have been completely or partially analyzed and some of them outperform the classical AS. Despite of many contributions, the analysis of this method remains incomplete.

We mention some of the developments related to the RAS method. The methods was introduced in [4]. The authors introduced the RAS as a cheaper and faster variants of the classical AS preconditioner for general sparse linear systems. The new method was shown to perform better that the AS according to the numerical studies presented there (see also [3]). The authors of [4] quoted that

…RAS was found accidentally. While working on a AS/GMRES algorithm in a Euler simulation, we removed part of the communication routine and surprisingly the “then AS” method converged faster both in terms of iteration counts and CPU time. We note that RAS is the default parallel preconditioner for nonsymmetric sparse linear systems in PETSc …

Many works have been devoted to RAS and therefore it would be difficult to present a complete review of them. Here we mention that in [7, 6] an algebraic convergence analysis is presented. In [2] the authors provide and extension of RAS using the so-called harmonic overlaps (RASHO). Both RAS and RASHO outperform their counterparts of the classical additive Schwarz variants. An almost optimal convergence theory is presented for the RASHO. In [5], it is shown that a matrix interpretation of RAS iteration can be related to the the continuous level of the underlying problem. The authors explain how this interpretation reveals why RAS converges faster than classical AS. Still, an explanation of the condition number of the RAS remains to be satisfactory. In [12], a by now classical book introducing domain decomposition methods, the authors comment

To our knowledge, a comprehensive theory of this algorithm is still missing. We note however that the restricted additive Schwarz preconditioner is the default parallel preconditioner for nonsymmetric systems in the PETSc library …and has been used for the solution of very large problems…

In this paper we re-visit the method proposed by Cai and Sarkis. We re-interpret the method as an iterative procedure where each iteration requires the solution of elliptic interface problems in each overlapping subdomain. The analysis of the method is presented in an abstract setting. First we write a Hilbert space framework for the analysis of the classical additive method. Then we generalize this Hilbert space framework and apply this extension to compute the condition number of several methods that use restrictions onto original subdomains in the construction (instead of restrictions to the overlapping subdomains). We present abstract results that may be useful to analyze non-symmetric domain decomposition method in general. We illustrate in particular how to use the results for a one level restricted additive method. Several other models and similar methods can be considered as well. For instance, restricted method for the elasticity equation, two-level domain decomposition method with classical or modern coarse spaces design, e.t.c.

The rest of the paper is organized as follows. In Section 2 we review the classical domain decomposition methods results in a simple Hilbert space framework. In Section 3 we recall the classical AS one level method. In Section 4 we present the abstract analysis of symmetric methods. We first revisit the analysis for symmetric methods using projections and angles between sub-spaces. We generalize this analysis to nonsymmetric methods. In particular we apply this analysis to a special family of nonsymmetric method. In Section 6 we define the restricted method that we analyze. In Section 7 we use the previously obtained results to write a condition number estimate of the restricted method defined before.

## 2 A Hilbert space framework

Let and be real Hilbert spaces with inner products and , respectively. The case of complex Hilbert spaces is similar. Consider to be a bounded operator with operator norm . In domain decomposition methods literature is referred to as a restriction operator. Introduce the transpose operator defined by

 a(RT,bϕ,v)=b(ϕ,Rv) for all v∈H. (1)

Despite of the fact that and operator norms depend on inner products and , our notation makes explicit only the dependence on . We use this convention also for operator norms.

Assume there is a (closed) subspace such that

 E=RT,b|G (2)

is easy to compute. The operator is known as an extension operator. Note also that and that we have with

 ET,b=ΠG,bR (3)

where is the orthogonal projection on using the inner product . We want to study the operator

 EET,b=RT,bΠG,bR. (4)

This operator is clearly symmetric and non-negative in the inner product. If we want to be non-singular and since and are orthogonal in , we need to be sure is 1-1 or, equivalently, is onto. A sufficient condition for the symmetric operator to be invertible is given by the following lemma known as stable decomposition lemma or Lion’s lemma in domain decomposition community. For the sake of completeness we show a detailed proof as it is usually presented in domain decomposition literature; see for instance [12, Chapter 2] or [9] and references therein. We note that we do not need to refer to the space

at this moment. Later we revise some of these inequalities in a more natural way to obtain a sharper estimate.

###### Lemma 1 (Lions Lemma)

Assume that there exists a bounded right inverse of . That is, there exists a bounded operator such that for all . Then, the mapping is non-singular. Moreover, we have

 ∥ˆE∥−2b∥v∥2a≤a(v,EET,bv)=∥ET,bv∥2b≤∥E∥2b∥v∥2a

for all .

Proof. Note that for we have,

 ∥v∥2a = a(v,v) = a(EˆEv,v) = b(ˆEv,ET,bv) ≤ b(ˆEv,ˆEv)1/2b(ET,bv,ET,b)1/2 ≤ ∥ˆEv∥ba(v,EET,bv)1/2 ≤ ∥ˆE∥b∥v∥aa(v,EET,bv)1/2.

Using this last inequality we obtain

 ∥ˆE∥−2b∥v∥2a≤a(v,EET,bv).

To obtain the upper bound we proceed as follows using properties of subordinated norm of operators,

 ∥EET,bv∥2a ≤ ∥E∥2b∥ET,bv∥2b = ∥E∥2bb(ET,bv,ET,bv) ≤ ∥E∥2ba(v,EET,bv) ≤ ∥E∥2b∥v∥a∥EET,bv∥a

and therefore . We also have,

 a(v,EET,bv) ≤ ∥v∥a∥EET,bv∥a≤∥E∥2b∥v∥2a.

This finishes the proof.

###### Remark 2

Note that what it is needed is the existence of operator such that is invertible. In this case we have that is an stable right inverse of .

If, in addition, the extension operator comes from a restriction operator , as in (2), we can state the following corollaries.

###### Corollary 3

Let be a restriction operator such that . Assume that there exits a bounded operator such that for all . Then, the mapping is non-singular with

 ∥JR∥−2b∥v∥2a≤a(v,EET,bv)=∥ET,bv∥2b≤∥R∥2b∥v∥2a

for all .

###### Corollary 4

Let be a restriction operator such that . Assume that there exits a bounded operator such that for all . Then, the mapping is non-singular with

 ∥ˆE∥−2b∥v∥2a≤a(v,RT,bRv)=∥Rv∥2b≤∥R∥2b∥v∥2a

for all .

Assume that is the solution of the following variational equation,

 a(u,v)=L(v) for all v∈H. (5)

Assuming that is easy to compute. We see that, for the solution , is possible to compute using this variational equation (without explicitly knowing or computing the function ). In fact, we have

 b(ET,bu,ϕ)=a(u,Eϕ)=L(Eϕ)=LEϕ for all ϕ∈H. (6)

This equation might be easier to solve numerically than the original problem. Therefore, we can alternatively compute the solution of (5) by iteratively solving the equation,

 EET,bu=˜L (7)

where . When implementing an iterative method, in each iteration we have to apply the operator

to a residual vector, say

. More precisely, we have to

1. Compute , this can be done by solving the equation

 b(x,ϕ)=a(r,Eϕ) for all ϕ∈H. (8)

In terms of the restriction operator we have .

2. Compute by applying the extension operator , that is whic is assumed possible and numerically efficient to compute.

The practicality of using this iteration depends on the possibility to inexpensively compute the right hand side and, of course, the condition number of . Once the right hand side is computed, then the performance of the iterative procedure depends on the condition number of the associated operator equation. If we use the spectral condition number of the operator , we see from Lemma 1 that

 κspectral(EET,b)≤∥ˆE∥2b∥E∥2b.

Then, the number of iterations for solving the equation (7) (up to a desired tolerance) will depend on . In this case, due to the symmetry and positivity of the operator we can use, for instance, a conjugate gradient algorithm to solve the equation above.

In general, the condition number of an operator is defined by

 κ(T)=∥T−1∥b∥T∥b.

Recall that we use the operator norm notation Here the norm where we make explicit the dependence on only.

## 3 Classical additive method for Laplace equation

In this section we use the Hilbert space framework above to review the analysis of the classical additive method. As usual, we consider a subdomain with a non-overlapping partition of the domain into subdomains . By enlarging these subdomains an specific width we obtain and overlapping decomposition . For more details see [12].

Let , and . In this case consider

 a(v,v)=∫D|∇v|2 for all v∈H.

Denoting by the elements of , we define,

 b({vℓ},{vℓ})=NS∑ℓ=1∫Oℓ|∇vℓ|2+∫∂Ωℓv2ℓ=NS∑ℓ=1bℓ(vℓ,vℓ)+b∂ℓ(vℓ,vℓ)

where we have put and .

###### Remark 5 (Norm boundary term)

The roll of is not essential and can be replaced by any other bilinear form that vanish for function on and makes a norm on .

Introduce also defined by . Equation (1) defining corresponds to

 ∫D∇RT,bvℓ∇z=b({vℓ},{zℓ|Oℓ}) for all z∈H.

Which implies that , defined as , is also given by

 E{vℓ}=NS∑ℓ=1Eℓvℓ

where is the extension by zero outside operator. To see this note that

 ∫D∇E{vℓ}∇z=NS∑ℓ=1∫Oℓ∇vℓ∇z for all z∈H.

We have that is given by

 b(ET,bu,{vℓ})=a(u,E{vℓ}).

Note that where solves the local equation

 ∫Oℓ∇ET,bℓu∇vℓ=∫D∇u∇Eℓvℓ for all vℓ∈H10(Oℓ).

Observe that can be obtained by solving a local problem. Then

In this case we denote . The existence of a right inverse can be stated as follows as it is common in domain decomposition literature.

Stable decomposition: There exists a constant such that for all there exist , such that and

 N∑ℓ=1bℓ(vℓ,vℓ)≤C2Ea(v,v). (9)

This clearly implies the existence of and . In fact, where the functions are the ones given by the stable decomposition assumption.

Strengthened Cauchy inequalities. There exits a matrix with and such that

 a(Eℓvℓ,Ekuk)≤μℓkbℓ(vℓ,vℓ)1/2bk(vk,vk)1/2.

By using bilinearity and vector Chauchy inequalities, this clearly implies that

 a(E{vℓ},E{vk}) = N∑ℓ,ka(Eℓu,Eku) ≤ ∑ℓ,kμℓ,kbℓ(vℓ,vℓ)1/2bk(vk,vk)1/2 ≤

where is the spectral radius of the matrix above. Then and using that and in Lemma 1 we have the following result.

###### Corollary 6

For all we have

where is the stable decomposition constant and is the spectral radius of the matrix above.

###### Remark 7

Let us consider the case of the one level additive method setting. More levels can be analyzed in a similar way. In the one level setting, with original domains of size and overlapping of size a usual bound for is given as follows by constructing a stable decomposition as follows; see [12]. Start by constructing cut of functions such that

 ηℓ(x)=1 for all x∈Dℓ, ηℓ(x)=0 % for all x∉Oℓ, |∇ηℓ(x)|≤Ccut1δ. (11)

Define the partition of unity function

 χℓ=ηℓη where η=NS∑ℓ=1ηℓ.

We see that

 χℓ(x)=1 for all x∈Dℓ∖∪ℓ′≠ℓOℓ′, ξℓ(x)=0 for all x∉Oℓ,

and therefore

 |χℓ(x)|≤Cpu1δ. (12)

We should have . Then define . Define

 Neigh(j)={j:Oj∩Oi≠∅} and ν=maxj#Neigh(j) (13)

We have (see [12])

 C2E⪯νCpu(1+1Hδ).

The norm of and are bounded by

 ∥E∥2≤ρ≤ν.

In case we want to approximate the solution of problem (5), we see that also solves the operator equation,

 EET,bu=(NS∑ℓ=1Pℓ)u=˜L.

Where is obtaining by -assembling the solutions of the local problems,

 bℓ(ET,bℓu,vℓ)=L(Eℓvℓ) for all vℓ∈H10(Oℓ).

The linear system in (7) is then well conditioned.

## 4 Nonsymmetric methods obtained by changing restrictions and the inner product

We use the Hilbert space framework introduced in Section 2. Recall that we have and Hilbert spaces with inner products and , respectively. We also used the bounded restriction operator for the construction.

To obtain nonsymmetric methods we additionally introduce a second bi-linear form defined on . Let us introduce a possibly different and bounded restriction operator and the transpose defined analogously to (1) by

 a(ST,cϕ,v)=c(ϕ,Sv) for all v∈H. (14)

Define as a second extension operator. As before assume that there is an stable left-inverse for , say such that for all (that is, is bounded in the and inner product norms). We can then conclude about and similar inequalities than the given before in the case is symmetric and positive definite. In particular is a bijective application from onto .

###### Corollary 8

Let be a restriction operator such that . Assume that there exits a bounded right inverse of , say . Then, the mapping is non-singular with

 ∥ˆF∥−2c∥v∥2a≤a(v,FFT,cv)=∥FT,cv∥2c≤∥F∥2c∥v∥2a

and

 ∥ˆF∥−2c∥v∥2a≤a(v,ST,cSv)=∥Sv∥2c≤∥S∥2c∥v∥2a

for all .

We want to study the nonsingularity of the operator . Note that

 FET,b=ST,cΠG,bR. (15)

As a particular case we can put . In this case, and . We can then obtain the operator

 FET,b=RT,cΠG,bR. (16)

See (4). This operator is also nonsymmetric for general bi-linear forms and . This is due to the fact that might not be symmetric in the bilinear form.

###### Remark 9 (Perturbation theory)

Note hat we can write

 EFT,b=EET,b+E(F−E)T,b=EET,b+J

where is a perturbation of of size . Several results can be pursue of the type: If is small, then the operator will be invertible and it is possible to estimate its condition number. Here we found these results are not practical for analyzing domain decomposition methods.

## 5 Condition number estimates using norms of projections

In this section we present a different analysis that may turn useful when estimating condition number of preconditioned operators (not-necessarily constructed by a domain decomposition design). We present a series of projection arguments in order to be able to study nonsymmetric methods. As presented earlier, the idea is to be able to estimate the condition number of an operator of the form where are different extension operators. In particular we are able to bound condition number for the family of nonsymmetric method presented in Section 4 where the extension operators are defined from restriction operator from to a bigger space . Before going to to nonsymetric method we revisit the condition number bound of symmetric methods.

### 5.1 Condition number of symmetric methods revisited

There is a simple way to interpret the non-singularity of obtained from the existence of , the right inverse of . We can construct solutions of the equation

 EET,bw=u (17)

which is equivalent to

 RT,bΠG,bRw=u (18)

as follows. The solution can be constructed by applying projections defined on . First recall that if we have , and , where is the identity operator on , we also have that . Therefore, has a left stable inverse that can be used in the analysis.

Let be given. Equation (17) is equivalent to

 a(EET,bw,v)=a(u,v) for all v∈H.

Note that so that we have

 b(ET,bw,ET,bv)=b(ˆEu,ET,bv) for all v∈H

and therefore we conclude that is the orthogonal projection of on . We can then construct as follows:

1. Define . Then we readily see that . By assumption we then have

 ∥y∥b≤∥ˆE∥b∥u∥a. (19)
2. Construct such that . In fact, we can use the orthogonal projection onto the subspace , which is denoted by . In this case we take . Note that . Therefore, this projection is along the subspace we have and therefore . In this case we have the obvious estimate

 ∥x∥b≤∥ΠR(ET,b),by∥b≤∥y∥b. (20)

See Figure 1 for an illustration.

3. Observe that and therefore for some . By applying we can make explicit to get . Then is the solution of (18). We obviously have the estimate

 ∥w∥b≤∥ˆET,b∥b∥x∥b=∥ˆE∥b∥x∥b. (21)

Combining (19), (20) and (21) give us

 ∥w∥a≤∥ˆE∥b∥x∥b≤∥ˆE∥∥y∥b≤∥ˆE∥2b∥u∥a.

We then obtain that is invertible and .

It is easy to see that finally given a bound for the condition number of as

 κ(EET,b)=∥EET,b∥b∥(EET,b)−1∥b≤∥ˆE∥2b∥E∥2b.

Note that this is not an spectral condition number but rather the condition number of the operator .

A first consequence of this analysis is that it is evident that the estimate in (20) is, in general, not sharp. It does not have into account the relative position of the subspace with respect to , which may be taken into account.

We need the following definitions and results; see [11, 1, 8]. Let and be subspaces of (or ). Introduce the minimal angle between subspaces and with respect to the inner product , as

 sin(θb(X,Y))=inf∥x∥=1,x∈Xdistb(x,Y)=inf∥x∥b=1,x∈X√1−∥ΠYx∥2b. (22)

Equivalently, we have where

 cos(θb(X,Y))=supx∈X,y∈Yb(x,y)∥x∥b∥y∥b=∥ΠXΠY∥b=∥ΠXΠY∥b. (23)

Still equivalent, we have,

 sin(θb(X,Y))=1∥Q(X,Y)∥b (24)

where is the (oblique) projection on and in the direction of .

Introduce the maximal angle between subspaces and , as

 sin(Θb(X,Y)) = sup∥x∥b=1,x∈Xdistb(x,Y) = sup∥x∥b=1,x∈X√1−∥ΠYx∥2b=∥ΠX−ΠY∥b.

Equivalently we have where

 cos(Θb(X,Y))=infx∈X,y∈Yb(x,y)∥x∥b∥y∥b. (25)

We also have,

 θb(X,Y)+Θb(X,Y⊥,b)=π2,

and

 sin(θb(X,Y))=cos(Θb(X,Y⊥,b)). (26)

Denote by the restriction of to . Then we can replace (20) by the sharper estimate

 ∥x∥b≤∥ΠR(ET,b)y∥b≤∥ΠR(ET)|R(ˆE)∥b∥y∥b. (27)

Observe that (see [11, 1, 8])

 ∥ΠR(ET,b)|R(ˆE)∥b = sup∥x∥b=1,x∈R(ˆE)∥ΠR(ET,b)x∥b (28) = sup∥x∥b=1,x∈R(ˆE)√1−∥ΠN(E)x∥2b (29) = sin(Θb(R(ˆE),N(E))) (30) = cos(θb(R(ˆE),R(ET,b))). (31)

See Figure 1 for an illustration. See [11, 1, 8] more details and related results on oblique projections.

###### Remark 10

It is clear that for the case of the classical method of Section 3 where the overlap parameter is small we have but in general (for wide for instance) we may have .

We then have the following somehow sharper result for the abstract case.

###### Lemma 11

Assume that there exists a bounded right inverse of . That is, there exists a bounded operator such that for all . Then, the mapping is non-singular. Moreover, we have

 ∥(EET,b)−1∥b≤cos(αE)∥ˆE∥2b

where is the minimal angle between subspaces and , that is,

 αE=θb(R(ˆE),R(ET,b)).

We then have the bound

 κ(EET,b)=∥EET,b∥b∥(EET,b)−1∥b≤cos(αE)∥ˆE∥2b∥E∥2b.

Proof. The estimate is obtained by combining (19), (21) and the bound (27).

###### Remark 12

The operator is obviously non-negative and therefore the condition number estimate in Lemma 11 we have useful bounds for converge of iterative method such as Krylov subspace methods (when is of finite dimension, for instance).

There is another interesting observation that is useful for the analysis and it is worth to state as a result before going to nonsymmetric methods.

###### Lemma 13

The operator is a projection on and along . Analogously, the operator is a projection on and along .

Using this lemma we can study the relative position of subspaces of interest. For instance we have,

 cos(Θb(R(ˆE),R(ET,b))=sin(θb(R(ˆE),N(E))=1∥QE∥b≥1∥ˆE∥b∥E∥b. (32)
###### Remark 14

Note that and are inverse to each other.

### 5.2 General nonsymmetric method analysis using projections

Let be a second extension operator. We want to study the operator . See Figure 2 for an illustration.

###### Theorem 15

Consider extensions operators and with stable right inverse and , respectively. Assume the boundedness of , the oblique projection onto and in the direction of . Then, the operator is invertible. Moreover,

 ∥(FET,b)−1∥b≤∥Q|R(ˆE)∥b∥ˆE∥b∥ˆF∥b.

Proof. As introduced before, we solve the equation

 FET,bw=u.

Let be given.

1. Define . Then we readily see that . By assumption we then have

 ∥y∥b≤∥ˆF∥b∥u∥a. (33)
2. Construct such that . Here we use the oblique projection . See Figure 2. In fact, . By definition of the projection we have so that . We have,

 ∥x∥b≤∥Q|R(ˆE)∥∥y∥b. (34)
3. Take such that . In fact, . This is the solution of the equation above since we have . We can bound

 ∥w∥b≤∥ˆE∥b∥x∥b. (35)

By combining the estimates in (33), (34) and (35) above we finish the proof.

We can now give a bound for the condition number of the operator .

###### Corollary 16

We have

 κ(FET,b)=∥FET,b∥b∥(FET,b)−1∥b≤∥ˆF∥b∥F∥b∥Q|R