Optimality properties of a square block matrix preconditioner with applications

Computers and Mathematics with Applications xxx (xxxx) xxx

Contents lists available at ScienceDirect

Computers and Mathematics with Applications journal homepage: www.elsevier.com/locate/camwa

Optimality properties of a square block matrix preconditioner with applications Owe Axelsson

∗

The Czech Academy of Sciences, Institute of Geonics, Ostrava, Czech Republic Department of Information Technology, Uppsala University, Uppsala, Sweden

article

info

a b s t r a c t A preconditioned square block matrix, called PRESB has previously been applied successfully and, for more standard type of problems, have been shown to have eigenvalue bounds in the interval (1/2, 1], which holds uniformly with respect to all parameters involved. Having such fixed bounds enables the use of an inner-product free acceleration method, such as the Chebyshev iterative method. Here it is shown that the method can be applied also for some more general problems, where the spectrum contains outlier eigenvalues larger than unity and that one can apply a polynomially modified version of the Chebyshev method for such problems. © 2019 Published by Elsevier Ltd.

Article history: Available online xxxx Keywords: Preconditioning Square block matrices Optimal control Subspace observation and control

1. Introduction Two-by-two block matrices with square blocks arise in many applications. An elementary example is when solving a complex valued system, (A + iB)(x + iy) = (f + ig), where A, B, f , g are real, in real arithmetics by rewriting it in two-by-two real block matrix form,

[ ] A

[

x A := y B

][ ] −B x A

y

[ ] =

f , g

(1)

see e.g. [1–3] for some early presentations of this method. We assume that A + B is nonsingular, but A and/or B can be singular. Another example arises in optimal control of solutions of partial differential equations, OPT-PDE, where in the simplest case a functional, 1 1 ∥u − ud ∥2 + β∥v∥2 2 2 is to be minimized under the constraint (again as a simple example only)

˜ , v) = J(u

−∆u + cu = v in Ω , where ∆ is the Laplacian operator and c ≥ 0 (here assumed to be constant), with boundary conditions, for simplicity chosen as homogeneous, u = 0 or ∂∂ un = 0 on ∂ Ω . Here Ω is a given bounded connected domain in Rd , d = 1, 2 or 3, u is the solution to be computed to become close to a given target, i.e. desired or observed solution ud , v is a control variable, ∗ Correspondence to: The Czech Academy of Sciences, Institute of Geonics, Ostrava, Czech Republic. E-mail address: [email protected]. https://doi.org/10.1016/j.camwa.2019.09.024 0898-1221/© 2019 Published by Elsevier Ltd.

Please cite this article as: O. Axelsson, Optimality properties of a square block matrix preconditioner with applications, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.09.024.

2

O. Axelsson / Computers and Mathematics with Applications xxx (xxxx) xxx

here distributed on the whole domain, chosen to achieve this goal and β > 0 is a regularization parameter which is needed to cope with the ill-conditioning, i.e. of an inverse operator type, of the problem. For simplicity we seek a solution in the Sobolev space, H01 (Ω ), i.e. H 1 (Ω ) with homogeneous Dirichlet conditions imposed, and let ∥ · ∥ be the L2 −norm. Such a problem has been considered in [4–6], see also [7] for more involved problems. ˜ , v ) is completed with a Lagrange multiplier, In order to solve the constrained optimization problem, the functional J(u that is, an adjoint variable, λ and a weakly regularized equation (here self-adjoint),

∫ Ω

λ(−∆u + cu − v ) =

∫ Ω

(∇λ · ∇ u + c λu − λv ).

The so completed Lagrange functional takes the form J(u, v, λ) =

1 2

∥ u − ud ∥ 2 +

∫ Ω

(∇λ · ∇ u + c λu − λv ) +

1 2

β∥v∥2 ,

where we seek the inf(u, v ) − sup(λ) solution. We restrict now the problem to a finite element space V0,h ⊂ H01 (Ω ), u ∈ Vh , v ∈ Vh ⊂ H 1 (Ω ). After application of the ∂J ∂J ∂J first order necessary (KKT) conditions, ∂ u = 0, ∂v = 0, ∂λ = 0, we get the block matrix algebraic system,

⎡

M ⎣0 ˜ K

0 βM −M

⎤[ ] [ ] ˜ u Mud K ⎦ 0 v = −M , 0 λ 0

(2)

where M is the mass matrix, ˜ K = K + cM and K is the stiffness matrix corresponding to the discretized Laplacian operator. As we shall see, there is a unique solution to this system. √ After elimination of v , the last equation scaled with β and introducing λˆ = − √1β λ, we get

] ][ ] [ u Mud −ˆ K , = 0 λˆ M √ where ˆ K = βˆ K , i.e. a system on the form (1). [

M ˆ K

(3)

One type of method that has been used to solve this type of equation is based on reduction of (3) to its Schur complement form, S u = (M + ˆ K M −1 ˆ K )u = Mud .

To avoid having to form the matrix explicitly, which is dense in general, one can use a factored preconditioner,

˜ = (M + ˆ S K )M −1 (M + ˆ K ), see e.g. [8] for an early presentation of such methods. The present analysis is done in a more compact and straightforward ˜ are transformed to I + ˜ way. Using a similarity transformation with M −1/2 , the matrices S and S K˜ K T respectively to T −1/2 −1/2 T ˜ ˜ ˜ ˜ ˜ (I + K )(I + K ), where K = M KM . Since (I − K )(I − K ) is symmetric and positive semi-definite, it follows that

˜ ˜−1 S are contained in the interval (1/2, 1], which shows that the condition K +˜ KT ≤ I +˜ K˜ K T . Hence the eigenvalues of S

number is bounded by 2. However this method requires that M is nonsingular since it needs a solution with M at each iteration step to evaluate the residual. In more general problems, M can be singular. The purpose of this paper is to present a preconditioning matrix for A and an efficient iterative acceleration method, the combination of which gives an iterative solution method of optimal order and a uniformly fast rate of convergence with respect to all parameters involved, including the mesh size parameter h. As we shall see, here appears no solution with matrix M. The method has been applied for several different types of optimal control problems for PDEs, see [2,4–7,9]. The present derivations are shorter and more general applicable. It gives sharp eigenvalue bounds which enables an efficient use of a Chebyshev iteration method. The paper is composed as follows. In next section we present the preconditioner and its spectral properties. In Section 3 we present a polynomially modified Chebyshev iteration method which enables use of the method also for problems where the eigenvalues are contained in two different intervals or there are outlier eigenvalues to the basic interval, which occurs for certain indefinte matrix problems. An example of such a problem is analysed in Section 4. This section contains also some numerical illustrations. 2. The preconditioner and the spectral properties of the preconditioned matrix 2.1. Spectral properties Consider first the system B

[ ] [ ξ A = η B

][ ] [ ] −B ξ 0 = , A + 2B η 2α Bη

∥ξ ∥ + ∥η∥ ̸= 0

(4)

where A and B are symmetric and positive semidefinite and A + B is nonsingular and α ̸ = 0. Please cite this article as: O. Axelsson, Optimality properties of a square block matrix preconditioner with applications, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.09.024.

O. Axelsson / Computers and Mathematics with Applications xxx (xxxx) xxx

Lemma 2.1.

3

Given (4), then

η = 2α G(I − G)η, where G = (A + B)−1 A. Proof. It holds Aξ = Bη

(5)

(A + B)(ξ + η) = 2α Bη.

(6)

and

Hence

ξ + η = 2α (A + B)−1 Bη = 2α (I − (A + B)−1 A)η and a multiplication with A and use of (5) show that (B + A)η = A(ξ + η) = 2α A(I − (A + B)−1 A)η, so

η = 2α G(I − G)η. □ Theorem 2.1. Let α = 1 in (4). It follows that

[

A B

] −B A

is nonsingular.

Proof. For α = 1 it follows from Lemma 2.1 that

η = 2G(I − G)η or ((I − G)2 + G2 )η = 0, which, since (I − G)2 and G2 are positive semidefinite and their kernels have the trivial intersection shows that η = 0 and from (6) hence also ξ = 0. □

[

]

A −B be a preconditioner to A. Then the eigenvalues λ of B−1 A satisfy λ = B A + 2B where µ is an eigenvalue of G = (A + B)−1 A.

Theorem 2.2. Let B =

1 (1 2

+ (1 − 2µ)2 ),

Proof. Consider

[ ] [ ] ξ ξ λB =A , η η

∥ξ ∥ + ∥η∥ ̸= 0,

or (1 − λ)B

[ ] [ ] ξ 0 = . η 2Bη

Here λ = 1 if and only if η ∈ N (B), any ξ . For η ∈ N (B)⊥ , i.e. λ ̸ = 1, with α = 1/(1 − λ) Lemma 2.1 shows that (1 − λ)η = 2G(I − G)η = 2

( =

1 2

( −2

1 2

(

1 2

)2 ) −G

( −

1 2

)) ( ( )) 1 1 −G + −G η 2

2

η,

or

λ=

1 2

(1 + (1 − 2µ)2 ). □

Corollary 2.1. If A and B are symmetric and positive semidefinite and N (A) ∩ N (B) = {∅}, then 1 2

≤

1 2

(1 + (1 − 2µ)2 )=λ ≤ 1,

where µ is an eigenvalue of (A + B)−1 A. If ˆ µ = µmax < 12 or 21 < ˆ µ = µmin ≤ 1 then

1 (1 2

+ (1 − 2ˆ µ)2 ) ≤ λ ≤ 1.

Please cite this article as: O. Axelsson, Optimality properties of a square block matrix preconditioner with applications, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.09.024.

4

O. Axelsson / Computers and Mathematics with Applications xxx (xxxx) xxx

Proof. This follows directly from Theorem 2.2. □ Remark 2.1. Clearly 0 ≤ µ ≤ 1 if yT Ay ≥ 0, that is, 12 ≤ λ ≤ 1. Further µ > 1 if ηT Aη < 0 but ηT (A + B)η > 0. Then 1 < λ ≤ (1 + (2µmax − 1)2 )/2. Finally, if ηT Aη < 0 and ηT (A + B)η < 0 then µ < 0 and 1 < λ ≤ (1 + (2µmax + 1)2 )/2. Consider now the case where A may be indefinite but A + B is still nonsingular. If A is indefinite then we assume that

α := max

ηT Bη <1 ηT (−A)η

for vectors η, where ηT Aη < 0. Note that if:

• ηT Aη ≥ 0, then 0 ≤ µ ≤ 1 1 • ηT Aη < 0, then µ = > 1. (Note the assumption ηT Bη 1− T η (−A)η

ηT Bη ηT (−A)η

< 1.)

It follows that (i) 1 ≥ λ ≥

⎧ ⎨1

2

+ 2( 12 − max µ)2 ,

⎩1,

if max µ <

1 2

,

if 1 ≥ max µ ≥

2

(ii) 1 ≤ λ ≤ λmax ≤

α = max η

1 2

+ (max µ − 12 )2 =

ηT Bη , ηT (−A)η

1 2

where ηT Aη ≥ 0 1 2

1 +α 2 + 2( 1−α − 12 )2 = 12 (1 + ( 11−α ) ), where

ηT Aη < 0.

It is seen that for vectors η where ηT Aη < 0, it is important that α is not too close to unity. As a conclusion, the eigenvalues λ of B−1 A are contained in two intervals [1/2, 1] and [1, λmax ]. 2.2. Spectral analyses for a subset observation domain The eigenvalue analysis for the model problem in Section 2.1 can be extended to problems where the observation domain is a nontrivial subset of Ω . Denote the corresponding mass matrix with M0 , which then replaces the 1,1 block in (2) and (3). To handle the spectral analysis for this case, we use an intermediate matrix

[ C=

˜ M K

−K T ˜ M

]

˜ = 1/2(M + M0 ). To find the eigenvalues of C −1˜ where M A we consider, µC

[ ] [ ] ξ ˜ ξ , =A η η

where

[ ˜= A

M0 K

−K T

]

M

[ =C+

] −E ξ . Eη

It follows that (1 − µ)(ξ T ξ + ηT η) = −ξ T E ξ + ηT E η where E = 1/2(M − M0 ) so

−γ ≤ 1 − µ ≤ γ where

γ = max

ηT (M − M0 )η 1 − λmin (M −1 M0 ) ηT E η = max = <1 ηT η ηT (M + M0 )η 1 + λmin (M −1 M0 )

and γ > 0 since λmax (M −1 M0 ) < 1, because Ω0 is a nontrivial subset of Ω . Hence 1 − γ ≤ µ ≤ 1 + γ . It follows that the condition number of B−1 C is bounded by (1 + γ )/(1 − γ ) and hence the eigenvalues of B−1 A = B−1 CC −1 A are bounded as 1 ≤ λ ≤ 1 + γ, 2(1 − γ ) which holds uniformly with respect to all parameters. Clearly γ → 0 as Ω0 approaches the whole domain Ω . Please cite this article as: O. Axelsson, Optimality properties of a square block matrix preconditioner with applications, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.09.024.

O. Axelsson / Computers and Mathematics with Applications xxx (xxxx) xxx

5

2.3. Implementation A factorization of

[

] −B , A + 2B

A B= B shows that

[ B=

][

I

0

−I

A+B

I 0

−B

][

I

A+B I

]

0 , I

that is,

[ B=

][

I

0 I

−I

I 0

0

][

A+B

][ −B A + B

I 0

I

0

][

0 I

I I

]

0 . I

Hence B−1 =

[

][

I

−I

0 I

(A + B)−1 0

][

0 I

I 0

][

B I

I 0

0 (A + B)−1

][

I I

0 I

]

so, besides a matrix vector multiplication with B and some vector additions, an action of B−1 involves two solutions of linear systems with matrix A + B. In many problems, A + B is better conditioned than A or B, which latter matrices can even be singular. These systems are normally solved by inner iterations which implies that a variable preconditioner modification of a conjugate gradient method see [10] or the flexible version of GMRES (see [11,12]) should be used. However, for fairly accurate inner solutions one can still use a Chebyshev outer iteration method, see Section 4 for numerical illustrations. As a good solution method for the inner systems with the matrix A + B one can often use an algebraic multigrid method, see e.g. [13]. Since that matrix is better conditioned when the regularization parameter β is small, the modified incomplete factorization method, [1,14] can work well also. 2.4. Estimates of number of iterations due to disjoint eigenvalues In many important applications and preconditioning methods there arise a more or less clustered spectrum, typically around unity after a proper scaling, and several larger eigenvalues, often spread over a long distance from unity. As has been shown e.g. in [1], if there is a cluster of eigenvalues bounded by lower and upper eigenvalue bounds a, b and p disjoint larger eigenvalues, it follows that for a conjugate gradient method applied for such a symmetric matrix, the number of iterations is bounded by k ≤ ln

2/

ε

ln σ −1 + p.

√

√

Here ε > 0 is the required relative residual tolerance, and σ = (1 − a/b)/(1 + a/b) see e.g. [1]. This bound indicates that at most p additional iterations are needed to cope with the larger eigenvalues. However, in practice, due to rounding errors this does not hold. As has been shown in [15] see also [1], the influence of rounding errors can be analysed by perturbing the spectrum by an amount resulting from loss of orthogonality of residual and search direction vectors. As has been shown e.g. in [1], if there are isolated eigenvalues smaller than a the resulting ∑qnumber of iterations increases in addition to number of disjoint eigenvalues typically with ln 2/ε replaced by ln 2/ε + i=1 ln b/λi , where {λi }, i = 1, . . . , q is the set of such small eigenvalues. Hence it is seen that it is better to use a preconditioning method which results in some larger additional eigenvalues then in a number of smaller eigenvalues. For instance, it is much better to use the modified incomplete factorization method [1,14,16], than the classic incomplete factorization method [17], which methods give eigenvalues larger than unity respectively small eigenvalues. Due to the sensitivity of conjugate gradient methods to rounding errors it can be of interest to examine the replacement of conjugate gradient methods as iterative acceleration methods with Chebyshev methods, which do not relay on use of orthogonalized search direction vectors. A further fact supporting this observation is that since there are no inner products needed for orthogonalization in this method, there is also no computational overhead for synchronization and global communication of inner products in Chebyshev iteration methods. Such communication overhead occurs when the method are implemented on distributed memory parallel computers and will be even more dominating on future generations of parallel computers. This is a consequence of that the steady exponential growth rate of computer chip density, accurately predicted by Moore’s law some 50 years ago, now have started to slow down. Since multicore and heterogeneous computer architectures are the presently only available hardware technologies, the only way to increase computational speed is to use more cores and parallel computers, but then the now nearly fixed communication overhead will stay as a fixed computational expense term. Therefore a performance improvement must come from use of software and algorithms, such as the Chebyshev type methods, that have less dominating communication overhead. Please cite this article as: O. Axelsson, Optimality properties of a square block matrix preconditioner with applications, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.09.024.

6

O. Axelsson / Computers and Mathematics with Applications xxx (xxxx) xxx

The problem with Chebyshev iterations is that in practice one needs accurate estimates of the eigenvalue bounds and even good estimates of some of the dominating eigenvalues in the spectrum. In this section we have shown that the PRESB preconditioning method is an excellent method to use for this purpose and in the next section we show also how the dominating eigenvalues can be estimated and implemented in a polynomially modified Chebyshev iteration method. 3. A modified Chebyshev iteration methods The Chebyshev iteration method, has been presented early in [18]. See also [1,12] for later presentations. It can be defined by the recursion 2

(

x(k+1) = αk−1 x(k) − x(k−1) − where x(−1) = 0, αk−1 = 1 −

(

b−a 2(b+a)

)

a+b

)2

r (k) + x(k−1) ,

k = 1, 2, . . . .

αk−1 , k = 1, 2, . . ., α0 = 0.

For problems with outlier eigenvalues on can first ‘kill’ them which can be done by first computing the eigenvector corresponding to the maximal eigenvalue, i.e. λmaxˆ x = B−1 Aˆ x. Let C = I − λ 1 AB−1 . Then CAˆ x = 0 and C b = C (b − Aˆ x), max so the corrected right hand side vector, 1

(

b˜ = I −

)

λmax

AB−1 b.

(7)

For the case of more outlier eigenvalues, say an eigenvalue λ2 , one can use an additional matrix factor to (7), 1

(

I−

λ2

)

B−1 A B−1 r .

Besides the computation of the dominating eigensolutions, the above method requires a separate solution with the preconditioner. One can also add some eigenvalues larger but close to unity to the interval [1/2, 1], used for the Chebyshev method, which will make the interval slightly larger but still preserve its rapid convergence property. If an eigenvalue is close to unity, it can be used to extend the interval [1/2, 1] to include this also. One can then apply the Chebyshev iteration method. However, due to numerical rounding errors, see [1,15] the dominating eigenvalues become awake again and one must repeat the elimination, i.e. ‘killing’, process. An example of application of this method follows in the next section. 4. An equation of motion As an example of algebraic equations with matrices in two-by-two block form with an indefinite matrix block, we consider the space discretized equation of motion in a direct frequency form, see [19,20]. Following [20] and some references therein, consider the n-dimensional degree of freedom (n-DOF) equation of motion: M q¨ + C q˙ + Kq = p, depending on time, where M , K are the mass and stiffness matrices, respectively, C is a viscous damping matrix and q, p are the configuration and dynamical force vectors, respectively. We assume that the Helmholtz type operator −ω2 M + K is nonsingular, that is, ω is either sufficiently small so that ω2 ηT M η < ηT K η or sufficiently large so that ω2 ηT M η > ηT K η for vectors where K − ω2 M is negative definite. We apply a driving circular force p(t) = f˜ eiωt with frequency ω, which leads to a steady-state solution, q(t) = q˜ (ω)eiωt where q˜ is the solution of

[−ω2 M + K + i(ωCV + CH )]˜q(ω) = f˜ where C = CV +

1

C , ω H

(8)

that is, includes both viscous and hysteretic damping matrices.

Since the coefficient matrix varies with the frequency, a direct solver must repeat the factorization for each frequency, which will be too costly so an iterative method must be used. Thereby, the use of an efficient preconditioner is crucial. To apply the PRESB preconditioner, we first rewrite (8) as (A + iB)(u + iv ) = f + ig , where A = K − ω2 M ,

B = ω CV + CH ,

q˜ = u + iv

and f˜ = f + ig ,

which is written in real valued form

[ A

u

−v

]

[ =

A −B

][

B A

u

−v

]

[ =

f

−g

]

.

Please cite this article as: O. Axelsson, Optimality properties of a square block matrix preconditioner with applications, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.09.024.

O. Axelsson / Computers and Mathematics with Applications xxx (xxxx) xxx

7

Fig. 1. Eigenvalue distributions for the preconditioned matrix; left: ω = π , right: ω = 4π . Note the logarithmic scale on the right graph.

Applying the eigenvalue analysis in Section 2 for the preconditioned matrix, B−1 A =

[ [

=

A −B

B

A + 2B

]

I 0

]−1 [

[

0 A − I −B

]

A −B

B A

]−1 [

B

A + 2B

0 0

]

0 , 2B

we find 1 2

≤λ≤

1 2

where α = maxη 1 2

≤

1 2

(

(

1+

1+α

)2 )

1−α

ηT (ωCV +CH )η , ηT (ω2 M −K )η

,

for vectors η where ηT (K − ω2 M)η < 0 (if any), and

(1 + (1 − 2 max µ)2 ) ≤ λ ≤ 1,

where max µ = maxη

ηT (ωCV +CH )η , ηT (−ω2 M +K )η

if ηT (K − ω2 M)η > 0 for all η.

We show now two typical eigenvalue distributions for the PRESB preconditioned matrix. Fig. 1 shows that the eigenvalues are mainly contained in the interval [1/2, 1], except for larger values of ω, where there appear some outlying eigenvalues larger than unity. Note the large number of eigenvalues equal to unity and the spread of the few outlier eigenvalues for ω = 4π . In order to give a detailed example of the appearance of eigenvalues and the application of the polynomially modified Chebyshev iteration method we follow [20] and consider the above example where we assume now specifically that the domain of definition Ω is the unit square, M = Im is the m’th order identity matrix, K is the model problem five-point difference matrix on a uniform grid with mesh size h = 1/(m + 1), CV = 10Im and CH = µK , with damping coefficient µ = 0.02. Making such specific choices enables the computation of the eigenvalues of A, B and A + B, which still gives insight into the behaviour of eigenvalues of B−1 A also for more general problems. Note then first that A and B have the same eigenvectors, sin kπ x sin lπ y, k, l = 1, 2, . . .. As is familiar and readily computed, the corresponding eigenvalues of K equal 4 kπ h lπ h (sin2 + sin2 ), k, l = 1, 2, . . . . h2 2 2 Since we are primarily interested in the problem for small values of h, we will use the limit values

λk,l (K ) =

lim λk,l (K ) = (k2 + l2 )π 2 ,

h→0

k, l = 1, 2, . . . .

We list the eigenvalues of the preconditioned matrix (B−1 A), for the two-by-two square-block matrix system. It follows then that the eigenvectors of A and B are equal and

λk,l (A)/π 2 = −˜ ω2 + k2 + l2 ,

k, l = 1, 2, . . .

Please cite this article as: O. Axelsson, Optimality properties of a square block matrix preconditioner with applications, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.09.024.

8

O. Axelsson / Computers and Mathematics with Applications xxx (xxxx) xxx Table 1 Eigenvalues λ(B−1 A) larger than unity.

ω/π

k

l

λ(A)/π

λ(B)/π

λ(A + B)/π

µmax

λ(B−1 A)

1 2 2 2

1 1 2 2

1 1 1 2

1 1 4

3.2 6.36 6.42 6.48

4.2 4.36 7.42 10.48

1/4.2 1.46 ≈ 0.9 ≈ 0.6

– – – 2.46

3 3 3

1 2 2

1 1 2

−7 −4 −1

9.52 9.58 9.64

2.52 5.58 8.64

3.78 1.72 1.12

22.0 3.5 1.27

4 4 4 4 4 4

1 2 2 3 3 3

1 1 2 1 2 3

−14 −11 −8 −6 −3

12.68 12.74 12.80 12.84 12.90

−1.32

−9.6

1.74 4.80 6.84 9.90

7.3 2.67 1.88 1.28

205 93 9.92 3.85 1.71 –

5 5 5 5 5 5 5 5 5 5

1 2 2 3 3 3 4 4 4 4

1 1 2 1 2 3 1 2 3 4

15.84 15.9 16.0 16.0 16.1 16.2 16.1 16.2 16.2 16.4

−7.16 −4.1 −1.0 −1.0

−2.21 −3.9 −16.0

λk,l (B)/π 2 =

10

π

−2

2

−23 −20 −17 −15 −12 −7 −8 −5 0 7

˜ ω + 0.02(k2 + l2 ) = 3.183˜ ω + 0.02(k2 + l2 ),

16.0 3.9 1.7 2.0 1.3 1.0 0.70

4.1 9.2 8.1 11.2 16.2 23.4

5.2 49.2 505 468 15.6 3.4 9.0 1.8 1 –

k, l = 1, 2, . . .

where ˜ ω = ω/π . Further, since A and B have the same eigenvectors, it follows that

λk,l (A + B) = λk,l (A) + λk,l (B). Hence the eigenvalues of (A + B)−1 B equal

µk,l =

λk,l (A) = λk,l (A + B) 1+

1 10

ω+0.02(k2 +l2 ) π ˜ −˜ ω2 +k2 +l2

,

k, l = 1, 2, . . .

and it follows from the above that the eigenvalues of B−1 A equal

λk,l (B−1 A) =

1 2

1

+ (1 − 2µk,l )2 , 2

k, l = 1, 2, . . .

In Table 1 we list its values larger than unity. The largest eigenvalue of B−1 A can be computed by evaluating the eigenvalue expressions given above for the smallest values of k, l. For more general problems this value can be used as an initial approximation in a power iteration method to compute λmax and the next following dominating eigenvalues. Given λmax we scale the preconditioned system, λ 1 B−1 (Auk − b), so that the maximal eigenvalue of the scaled matrix max becomes equal to unity. Let ˆ λ be the scaled eigenvalues. A polynomially modified Chebyshev iteration method can then be defined by use of the polynomial q2 (ˆ λ) = ˆ λ(2 − ˆ λ). Here q′2 (1) = 0. The scaled eigenvalues, corresponding to the initial interval

ˆ λ2 =

1

λmax

1 2

≤ λ ≤ 1, become now small, ˆ λ1 =

1 , 2λmax

. Hence

q2 (ˆ λ2 ) q2 (ˆ λ1 )

2−

=2

2−

1

λmax 1 2λmax

< 2.

Further let ˆ λ3 = 0.3 (This can be seen as a reasonable approximation of the smallest eigenvalues in the set of large eigenvalues.) Then q2 (1) q2 (ˆ λ3 )

=

1 0.3 · 1.7

< 2.

Please cite this article as: O. Axelsson, Optimality properties of a square block matrix preconditioner with applications, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.09.024.

O. Axelsson / Computers and Mathematics with Applications xxx (xxxx) xxx

9

It follows that the polynomially modified Chebyshev polynomial, Tn (q2 (ˆ λ)) can be defined by the eigenvalue ratio a/b = 1/2 and that it converges well for both sets of unscaled eigenvalues, 1/2 ≤ λ ≤ 1, and the interval min{λ of eigenvalues larger then unity} ≤ λ ≤ λmax . It has been shown that for the direct frequency form of the equation of motion, the eigenvalues of the PRESB preconditioned matrix contain outlier eigenvalues above unity. These can be handled in two steps. The very large values are eliminated as shown in Section 3. The remaining eigenvalues are then contained in two intervals and a simple quadratic polynomial transforms them to a single interval, for which the standard but polynomially modified Chebyshev method can be applied. Hence the Chebyshev iteration method converges as fast as for the case of a single interval, (1/2, 1], and the increased computational cost comes only from the increased amount of matrix vector multiplications at each iteration, which is nearly doubled. Acknowledgements The computation of the eigenvalue distribution in Fig. 1 has been done by Zhao-Zheng Liang, which is much appreciated. Variable comments by two reviewers are also acknowledged. This work was supported by The Ministry of Education, Youth and Sports of the Czech Republic from the National Programme of Sustainability (NPU II), project ‘‘IT4Innovations excellence in science - LQ1602’’. References [1] O. Axelsson, Iterative Solution Methods, Cambridge University Press, Cambridge, 1994. [2] O. Axelsson, M. Neytcheva, Algebraic multilevel iteration method for Stieltjes matrices, Numer. Linear Algebra Appl. 1 (1994) 213–236, http://dx.doi.org/10.1002/nla.1680010302. [3] O. Axelsson, A. Kucherov, Real valued iterative methods for solving complex symmetric linear systems, Numer. Linear Algebra Appl. 7 (2000) 197–218, http://dx.doi.org/10.1002/1099-1506(200005)7:4<197::AID-NLA194>3.0.CO;2-S. [4] O. Axelsson, S. Farouq, M. Neytcheva, Comparison of preconditioned Krylov subspace iteration methods for PDE-constrained optimization problems: Stokes control, Numer. Algorithms 74 (2017) 19–37, http://dx.doi.org/10.1007/s11075-016-0136-5. [5] O. Axelsson, S. Farouq, M. Neytcheva, Comparison of preconditioned Krylov subspace iteration methods for PDE-constrained optimization problems, Poisson and convection–diffusion control, Numer. Algorithms 73 (2016) 631–663, http://dx.doi.org/10.1007/s11075-016-0111-1. [6] O. Axelsson, S. Farouq, M. Neytcheva, A preconditioner for optimal control problems, constrained by Stokes equation with a time-harmonic control, J. Comput. Appl. Math. 310 (2017) 5–18, http://dx.doi.org/10.1016/j.cam.2016.05.029. [7] O. Axelsson, D. Lukáš, Preconditioning methods for eddy-current optimally controlled time-harmonic electromagnetic problems, J. Numer. Math. 27 (2019) 1–21, http://dx.doi.org/10.1515/jnma-2017-0064. [8] J.W. Pearson, A.J. Wathen, A new approximation of the Schur complement in preconditioners for PDE constrained optimization, Numer. Linear Algebra Appl. 19 (2012) 816–829, http://dx.doi.org/10.1002/nla.814. [9] Z.-Z. Liang, O. Axelsson, M. Neytcheva, A robust structured preconditioner for time-harmonic parabolic optimal control problems, Numer. Algorithms 79 (2018) 575–596, http://dx.doi.org/10.1007/s11075-017-0451-5. [10] O. Axelsson, P.S. Vassilevski, Variable-step multilevel preconditioning methods, I: Self adjoint and positive definite elliptic problems, Numer. Linear Algebra Appl. 1 (1994) 75–101, http://dx.doi.org/10.1002/nla.1680010108. [11] Y. Saad, M.H. Schultz, GMRES: a generalized minimal residual algorithm for solving nonsymmetric linear systems, SIAM J. Sci. Stat. Comput. 7 (1986) 856–869, http://dx.doi.org/10.1137/0907058. [12] Y. Saad, Iterative Methods for Sparse Linear Systems, SIAM, Philadelphia, PA, 2003. [13] Y. Notay, An aggregation-based algebraic multigrid method, ETNA 37 (2010) 123–146. [14] I. Gustafsson, A class of first-order factorization methods, BIT 18 (1978) 142–156, http://dx.doi.org/10.1007/BF01931691. [15] A. Greenbaum, V. Pták, Z. Strakoš, Any nonincreasing convergence curve is possible for GMRES, SIAM J. Matrix Anal. Appl. 17 (1996) 465–469, http://dx.doi.org/10.1137/S0895479894275030. [16] O. Axelsson, Solution of linear systems of equations: Iterative methods, in: V.A. Barker (Ed.), Sparse Matrix Techniques, in: Lecture Notes in Mathematics, vol. 572, Springer, Berlin, Heidelberg, 1977, http://dx.doi.org/10.1007/BFb0116614. [17] J.A. Meijerink, H. van der Vorst, An iterative solution method for linear systems of which the coefficient matrix is a symmetric M-matrix, Math. Comp. 31 (1977) 148–162, http://dx.doi.org/10.1090/S0025-5718-1977-0438681-4. [18] R.S. Varga, Matrix Iterative Analysis, Prentice-Hall, Inc, Englewood Cliffs, N. J, 1962. [19] A. Feriani, F. Perotti, The formation of viscous damping matrices for the dynamic analysis of MDOF systems, Earth Eng. Struct. Dyn. 25 (1996) 689–709, http://dx.doi.org/10.1002/(SICI)1096-9845(199607)25:7<689::AID-EQE575>3.0.CO;2-L. [20] A. Feriani, F. Perotti, V. Simoncini, Iterative system solvers for the frequency analysis of linear mechanical systems, Comput. Methods Appl. Mech. Engrg. 190 (2000) 1719–1739, http://dx.doi.org/10.1016/S0045-7825(00)00187-0.

Please cite this article as: O. Axelsson, Optimality properties of a square block matrix preconditioner with applications, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.09.024.