This paper describes a multilevel preconditioning technique for solving complex symmetric sparse linear systems. The coefficient matrix is first decoupled by domain decomposition and then an approximate inverse of the original matrix is computed level by level. This approximate inverse is based on low rank approximations of the local Schur complements. For this, a symmetric singular value decomposition of a complex symmetric matix is used. Using the example of Maxwell's equations the generality of the approach is demonstrated

Schlundt, Rainer

English

Publications Server of the Weierstrass Institute for Applied Analysis and Stochastics

Weierstraß-Institutfür Angewandte Analysis und StochastikLeibniz-Institut im Forschungsverbund Berlin e. V.Preprint ISSN 2198-5855A multilevel Schur complement preconditioner for complexsymmetric matricesRainer Schlundtsubmitted: November 29, 2017Weierstrass InstituteMohrenstr. 3910117 BerlinGermanyE-Mail: rainer.schlundt@wias-berlin.deNo. 2452Berlin 20172010 Mathematics Subject Classification. 65F08, 65F15, 65N22, 65Y05.Key words and phrases. Complex symmetric sparse linear system, Schur complement, multilevel preconditioner, domaindecomposition, low rank approximation.Edited byWeierstraß-Institut für Angewandte Analysis und Stochastik (WIAS)Leibniz-Institut im Forschungsverbund Berlin e. V.Mohrenstraße 3910117 BerlinGermanyFax: +49 30 20372-303E-Mail: preprint@wias-berlin.deWorld Wide Web: http://www.wias-berlin.de/A multilevel Schur complement preconditioner for complexsymmetric matricesRainer SchlundtAbstractThis paper describes a multilevel preconditioning technique for solving complex symmetricsparse linear systems. The coefficient matrix is first decoupled by domain decomposition andthen an approximate inverse of the original matrix is computed level by level. This approximateinverse is based on low rank approximations of the local Schur complements. For this, a sym-metric singular value decomposition of a complex symmetric matix is used. Using the example ofMaxwell’s equations the generality of the approach is demonstrated.1 IntroductionWe consider iterative methods for solving large sparse systemsAx = b, (1)where A ∈ Cn×n, A = AT , A 6= AH , b ∈ Cn, and x ∈ Cn. Krylov subspace methods combinedwith a preconditioner solve the above system (1). For example, left preconditioning consists of modify-ing the original system into the system M−1Ax = M−1b. The preconditioner M is an approximationto A. The solve of the preconditioned system is relatively inexpensive.The domain decomposition (DD) approach decouples the original matrix A. We do not form the globalSchur complement system and do not solve it exactly. Let A be partitioned in 2× 2 block form asA =(B EET C), (2)where B ∈ Cm×m, C ∈ Cs×s, E ∈ Cm×s, and n = m+ s. We will receive the following basic blockfactorization of (2)A =(B EET C)=(I 0ETB−1 I)(B 00 S)(I B−1E0 I), (3)where S ∈ Cs×s, S = C −ETB−1E, is the Schur complement. UsingA−1 =(I −B−1E0 I)(B−1 00 S−1)(I 0−ETB−1 I), (4)the original system (1) can be easily solved if S−1 is available. The goal is to show that S−1 ≈C−1 + LRA = S̃−1, where LRA stands for low rank approximation matrix. The preconditioner Mthen has the following formM =(I 0ETB−1 I)(B 00 S̃)(I B−1E0 I). (5)DOI 10.20347/WIAS.PREPRINT.2452 Berlin 2017R. Schlundt 2We can writeS = C −ETB−1E = C1/2(I − C−1/2ETB−1EC−1/2)C1/2 = C1/2(I −G)C1/2 (6)andS−1 = C−1/2(I −G)−1C−1/2 = C−1 + C−1/2G(I −G)−1C−1/2 . (7)The symmetric matrix G ∈ Cs×s has a symmetric singular value decomposition (SSVD)G = C−1/2ETB−1EC−1/2 = WΣW T , (8)where W is a unitary matrix and Σ = diag(σ1, . . . , σs) with nonnegative σi (cf. [1]). Then LRA is anapproximation of C−1/2G(I−G)−1C−1/2. Detailed information can be found in Section 3. In Section2, we introduce the domain decomposition framework. Numerical experiments of a model problem arepresented in Section 4.Implementation details for a real symmetric matrix A are described in [4, 5].2 Domain decomposition frameworkAn interesting class of domain decomposition methods is the hierarchical interface decomposition(HID) ordering (cf. [2]). An HID ordering can be obtained from a standard graph partitioning (cf.METIS [3]). The reordered matrix has the following multilevel recursive form :Aj = PjCj−1PTj =(Bj EjETj Cj)and C0 ≡ A for j = 1, . . . , L . (9)Pj is a permutation matrix and L the number of levels. Each block Bj in Aj has a block-diagonalstructure resulting from this HID ordering. Analogous to (2), let Aj be partitioned at level j in blockform asAj =(Bj EjETj Cj)=Bj1 Ej1. . ....Bjp EjpETj1 . . . ETjp Cj, (10)where Bj ∈ Cmj×mj is a block-diagonal matrix, Bj = diag(Bj1 , . . . , Bjp), Cj ∈ Csj×sj , Ej ∈Cmj×sj , ETj = (ETj1, . . . , ETjp), Bji ∈ Cmji×mji , Eji ∈ Cmji×sj , 1 ≤ i ≤ p, nj = mj + sj , andmj = mj1 + · · ·+mjp . The following diagram illustrates these dependencies:(B1 E1ET1 C1)P2−→(B2 E2ET2 C2)P3−→ · · ·PL−1−−−→(BL−1 EL−1ETL−1 CL−1)PL−→(BL ELETL CL).Algorithm 1 provides a detailed illustration of the multilevel HID ordering scheme.3 Multilevel preconditioning techniqueAnalogous to (3), at each level j, the factorization of Aj is determined byAj =(Bj EjETj Cj)=(I 0ETj B−1j I)(Bj 00 Sj)(I B−1j Ej0 I), (11)DOI 10.20347/WIAS.PREPRINT.2452 Berlin 2017Schur complement preconditioner 3Algorithm 1 Basic pseudocode of multilevel HID orderingInput: Matrix A and maximum level L.Output: Matrices Bj , Cj , and Ej for 1 ≤ j ≤ L.1: procedure MHID2: Set C0 = A and s0 = n.3: for j = 1, . . . , L do4: Compute permutation matrix Pj and nj = sj−1.5: Re-sort matrix Cj−1 by Aj = PjCj−1PTj =(Bj EjETj Cj).6: nj = mj + sj . ⊲ mj = mj1 + · · ·+mjp7: Compute Bj , Bj = diag(Bj1 , . . . , Bjp). ⊲ Bj ∈ Cmj×mj , Bji ∈ Cmji×mji8: Compute Cj . ⊲ Cj ∈ Csj×sj9: Compute Ej , ETj = (ETj1, . . . , ETjp). ⊲ Ej ∈ Cmj×sj , Eji ∈ Cmji×sj10: end for11: end procedurewhere Sj = Cj −ETj B−1j Ej is the Schur complement at level j. ThusA−1j =(I −B−1j Ej0 I)(B−1j 00 S−1j)(I 0−ETj B−1j I)(12)is the inverse of Aj . Analogous to (7), S−1j can be approximated by C−1j plus an approximation ofC−1/2j Gj(I −Gj)−1C−1/2j . The preconditioner Mj then has the following formMj =(I 0ETj B−1j I)(Bj 00 S̃j)(I B−1j Ej0 I)(13)andC−1j = PTj+1M−1j+1Pj+1 . (14)At each level j, the symmetric matrix Gj ∈ Csj×sj has a singular value decomposition (SVD)Gj = C−1/2j ETj B−1j EjC−1/2j = UjΣjVHj , (15)where Uj and Vj are unitary matrices and Σj = diag(σj1, . . . , σjsj ) the singular values with realnonnegative σji . For the matrix Gj there exists a unitary matrix Wj such thatGj = C−1/2j ETj B−1j EjC−1/2j = WjΣjWTj (16)is an SSVD. An SSVD of a symmetric matrix can be determined from its SVD. Therefore we have tomodify the singular vectors corresponding to nonzero singular values (cf. [1]).Using the Sherman-Morrison-Woodbury (SMW) formula for (I − Gj)−1 the matrix Gj(I − Gj)−1results inGj(I −Gj)−1= WjΣjWTj(I +Wj(I − ΣjWTj Wj)−1ΣjWTj)= WjΣjWTj(W−Tj Σ−1j(I − ΣjWTj Wj)W−1j + I)·Wj(I − ΣjWTj Wj)−1ΣjWTj= Wj(I − ΣjWTj Wj)−1ΣjWTj(17)DOI 10.20347/WIAS.PREPRINT.2452 Berlin 2017R. Schlundt 4and consequentlyC−1/2j Gj(I −Gj)−1C−1/2j = C−1/2j Wj(I − ΣjWTj Wj)−1ΣjWTj C−1/2j= Zj(I − ΣjZTj CjZj)−1ΣjZTj(18)with Zj = C−1/2j Wj .The matrix C−1j ETj B−1j EjC−1j has the same singular values as C−1/2j ETj B−1j EjC−1/2j and thecolumns of C−1/2j Wj in (18) are the singular vectors of C−1j ETj B−1j EjC−1j due to the followingrelationGj = WjΣjWTj ⇐⇒ C−1/2j GjC−1/2j = ZjΣjZTj . (19)Thus, the computation of a low rank approximation to S−1j −C−1j (cf. (7), (18)) can be obtained by thefollowing SSVD problemC−1j ETj B−1j EjC−1j = ZjΣjZTj . (20)Let Σ̃j be the k largest singular values of (20) and Z̃j the corresponding singular vectors thenS−1j − C−1j ≈ Z̃j(I − Σ̃jZ̃Tj CjZ̃j)−1Σ̃jZ̃Tj (21)is a low rank approximation. Algorithm 2 provides a construction scheme of multilevel Schur comple-ment preconditioners based on low rank approximations.Algorithm 2 Construction of multilevel Schur complement preconditionresInput: Maximum level L, rank k, and matrices Bj , Cj , and Ej for 1 ≤ j ≤ L.Output: Matrices Z̃j and Σ̃j for 1 ≤ j ≤ L.1: procedure MSCP2: Approximate C−1L . ⊲ CL → C−1L3: for j = L, . . . , 1 do4: Approximate B−1j = diag(B−1j1, . . . , B−1jp ). ⊲ Bj → B−1j5: Compute the k largest singular values and the corresponding singular vectors fromC−1j ETj B−1j EjC−1j = ZjΣjZTj .⊲ Call Algorithm 3 to apply C−1j6: Set Σ̃j ← Σj and Z̃j ← Zj .7: end for8: end procedureUsing Eq. (14) we get the product of C−1j with a vector vj in the following way:yj = C−1j vj = PTj+1M−1j+1Pj+1vj = PTj+1M−1j+1uj+1 = PTj+1xj+1 .A recursive scheme for computing the product xj = M−1j uj is described in Algorithm 3. The solutionsassociated with Bj can be performed independently for each diagonal block in Bj .Finally, the preconditioned systemM−1Ax = M−1b with M−1 = C−10 = PT1 M−11 P1 (22)is to be solved.DOI 10.20347/WIAS.PREPRINT.2452 Berlin 2017Schur complement preconditioner 5Algorithm 3 A recursive formula for the approximation of xj = M−1j ujInput: Matrix Aj , vector uj at level j, and maximum level L.Output: Vector xj .1: procedure RSC(j,Aj ,uj)2: Split uj = (uj,1, uj,2)T . ⊲ uj = Pjuj−1,23: Compute z1 = B−1j uj,1.4: Compute z2 = uj,2 −ETj z1.5: if j = L then6: Compute xj,2 = C−1j z2.7: else8: Compute v = Pj+1z2.9: Compute w =RSC(j + 1,Aj+1,v).10: Compute xj,2 = PTj+1w.11: end if12: Compute xj,2 = xj,2 + Z̃j(I − Σ̃jZ̃Tj CjZ̃j)−1Σ̃jZ̃Tj z2. ⊲ c.f. Eq. (21)13: Compute xj,1 = z1 −B−1j Ejxj,2.14: return xj = (xj,1, xj,2)T . ⊲ xj−1,2 = PTj xj15: end procedure4 Numerical experimentsUsing the example of Maxwell’s equations we demonstrate the generality of the approach. We obtainin vector notation the following equations in integral form:˛P~E · d~l = −∂∂t¨A~B · d ~A˛P~H · d~l =∂∂t¨A~D · d ~A+¨A~J · d ~A‹S~B · d~S = 0‹S~D · d~S =˚Vq dV .(23)The constitutive relations belonging to them are~D = ε ~E , ~B = µ ~H , ~J = κ ~E . (24)Here, A is a surface with boundary curve P , V is a volume bounded by a surface S, and q is thevolume charge density. An orthogonal dual mesh is used to discretize the Maxwell’s equations usingthe Finite Integration Technique (FIT, [9, 10, 6]). The electric and magnetic voltages and fluxes overelemtary objects are defined as state variables in the following way:ei =ˆLi~E · d~l hj =ˆL̃j~H · d~l i = 1, . . . , nedi =¨Ãi~D · ~n d ~A bj =¨Aj~B · ~n d ~A j = 1, . . . , nfji =¨Ãi~J · ~n d ~A qk =˚Ṽkq dV k = 1, . . . , np .DOI 10.20347/WIAS.PREPRINT.2452 Berlin 2017R. Schlundt 6where ~n is the outward-pointing normal of the faces Aj and Ãi, respectively. If all field quantities varysinusoidally with time, the coefficient matrices of the corresponding linear systems of equation arecomplex, symmetric, and indefinite. Using Krylov subspace methods, (22) can be solved iteratively(cf. [7, 8]).We consider different dimensions of the coefficient matrices of the corresponding systems of linearequations, in fact n = 16632, n = 40824, and n = 472416. At each level j, j = 1, 2, . . . , thematrix A is partitioned into p, p ∈ {0, 2, 3, 5, 10, 15}, non-overlapping subsets Bji, i = 1, . . . , p. Ateach p, p ∈ {2, 3, 5, 10, 15}, we compute the k, k ∈ {5, 10, 15, 20}, largest singular values and thecorresponding singular vectors to obtain a low rank approximation. For p = 0 the solution is computedby [7, 8]. The same applies for the computation of the solution with the coefficient matrices Bji . FromTables 1-3, we find the numbers of iterations for the different dimensions. From Tables 4-6, we findthe corresponding dimensions of the matrices Bji and Cj for j = 1, 2, . . . and i = 1, . . . , p of theconsidered dimensions n = 16632, n = 40824, and n = 472416, respectively. It can be seen thatthe best results have been achieved for p = 2 and k ∈ {5, 10, 15, 20}. For small dimensions, alsop = 3 is useful. For n = 16632 and n = 40824, respectively, it is useful to choose k ∈ {5, 10}.For greater dimensions k ∈ {15, 20} is used. Experimental results indicate that this preconditionerbased on Schur complement approach is robust and can achieve savings in the iteration phase.Table 1: The number of iterations for n = 16632.number of k largest singular valuessubsets (p) 5 10 15 200 1572 128 124 126 1313 142 147 142 1435 146 143 143 14310 182 187 187 18015 211 208 217 214Table 2: The number of iterations for n = 40824.number of k largest singular valuessubsets (p) 5 10 15 200 5102 172 168 173 1673 182 173 176 1815 294 283 285 27510 347 341 330 29815 382 366 362 368DOI 10.20347/WIAS.PREPRINT.2452 Berlin 2017Schur complement preconditioner 7Table 3: The number of iterations for n = 472416.number of k largest singular valuessubsets (p) 5 10 15 200 9182 344 310 262 2703 537 532 550 5305 596 614 616 61710 716 758 719 76315 827 820 831 991Table 4: The dimensions of the matrices Bji and Cj for n = 16632.p level j dim(Bji) dim(Cj)21 7775 7775 10822 489 490 1033 41 41 2131 5121 4798 4915 17982 528 467 543 2603 82 76 68 3451 2913 3326 2793 2700 2770 21302 358 387 377 378 278 3523 57 47 63 62 47 76101 1663 1663 1336 1266 1094 40171161 1126 1238 1145 9232 310 280 310 247 258 1201333 300 202 284 2923 88 94 91 98 97 235107 95 110 90 96151 1108 1109 1109 479 787 5026804 605 836 636 668835 577 634 665 7542 252 254 214 205 284 1841148 163 209 177 218233 221 263 176 1683 94 87 92 101 81 53099 88 83 86 88108 71 88 55 904 29 24 14 17 22 19427 16 16 13 3120 25 27 27 28DOI 10.20347/WIAS.PREPRINT.2452 Berlin 2017R. Schlundt 8Table 5: The dimensions of the matrices Bji and Cj for n = 40824.p level j dim(Bji) dim(Cj)21 19324 19351 21492 967 973 2093 97 99 1331 12326 12174 12098 42262 1319 1279 1282 3463 103 109 103 3151 7018 7251 7043 5718 7133 66612 1213 1182 1163 1096 1146 8613 165 160 156 141 152 87101 3235 3201 3079 2481 3025 87893483 4082 3250 3045 31542 731 659 657 671 668 2052716 726 693 644 5723 184 166 191 156 186 301190 169 173 168 1684 28 24 26 24 24 2429 31 31 30 30151 2435 2430 2722 1894 2017 104041691 1979 1893 1957 16912048 2014 2069 1642 18482 570 537 456 533 408 2828512 530 593 469 361520 514 518 518 5373 161 154 159 132 158 586154 141 150 148 142141 139 142 166 1554 34 31 29 32 35 8238 35 36 33 3540 28 27 37 34DOI 10.20347/WIAS.PREPRINT.2452 Berlin 2017Schur complement preconditioner 9Table 6: The dimensions of the matrices Bji and Cj for n = 472416.p level j dim(Bji) dim(Cj)21 230148 230146 121222 5930 5935 25731 151607 150527 150873 194092 6351 6216 6239 6033 197 199 193 1451 89225 89134 89376 88286 85534 308612 5945 6049 6064 5930 5931 9423 188 187 187 189 183 8101 41695 43598 41198 43099 40019 5216741313 42221 42149 41569 433882 4818 4799 4870 4716 4929 36614910 4880 4709 4934 49413 361 359 360 356 361 65363 360 361 359 356151 27601 25284 26192 25952 26280 6894328966 27093 28287 28493 2650126212 25965 26983 26682 269822 4280 4186 4307 4297 4299 60204146 3961 4279 4262 40874239 4213 3869 4320 41783 391 376 391 390 392 209381 391 387 387 382388 396 387 392 3804 12 14 14 13 10 813 13 14 13 1414 15 13 13 16DOI 10.20347/WIAS.PREPRINT.2452 Berlin 2017R. Schlundt 105 ConclusionsThis paper presents a preconditioning method based on a Schur complement approach with low rankapproximations for solving complex symmetric sparse linear systems. This method can be both recur-sive and non-recursive. It tries to approximate the inverse of the Schur complement by exploiting lowrank approximations. For this, a hierarchical graph decomposition reorders the matrix into a multilevelblock form. On the negative side, building this preconditioner can be time consuming. A solve with thematrix Bj amounts to p local and independent solves with the matrices Bji, i = 1, . . . , p. These canbe carried out by a preconditioned Krylov subspace iteration. A big part of the computations to build apreconditioner based on Schur complement approach is attractive for massively parallel machines.References[1] Angelika Bunse-Gerstner and William Gragg. Singular value decomposition of complex symmet-ric matrices. Journal of Computational and Applied Mathematics, 21:41–54, 1988.[2] Pascal Henon and Yousef Saad. A parallel multistage ILU factorization based on a hierrarchicalgraph decomposition. SIAM J. Sci. Comput., 28(6):2266–2293, 2006.[3] G. Karypis and V. Kumar. A fast and high quality multilevel scheme for partitioning irregulargraphs. SIAM J. Sci. Comput., 20:359–392, 1998.[4] Ruipeng Li and Yousef Saad. Low-rank correction methods for algebraic domain decompositionpreconditioners. Technical Report ys-2014-5, Dept. Computer Science and Engineering, Univer-sity of Minnesota, Minneapolis, MN, 2014.[5] Ruipeng Li, Yuanzhe Xi, and Yousef Saad. Schur complement based domain decompositionpreconditioners with low-rank corrections. Technical Report ys-2014-3, Dept. Computer Scienceand Engineering, University of Minnesota, Minneapolis, MN, 2014.[6] Rainer Schlundt. Regular triangulation and power diagrams for Maxwell’s equations. WIASPreprint No. 2017, Weierstraß-Institut für Angewandte Analysis und Stochastik, 2014.[7] Rainer Schlundt, Franz-Josef Schückle, and Wolfgang Heinrich. Shifted linear systems in electro-magnetics. Part I: Systems with identical right-hand sides. WIAS Preprint No. 1420, Weierstraß-Institut für Angewandte Analysis und Stochastik, 2009.[8] Rainer Schlundt, Franz-Josef Schückle, and Wolfgang Heinrich. Shifted linear systems in electro-magnetics. Part II: Systems with multiple right-hand sides. WIAS Preprint No. 1646, Weierstraß-Institut für Angewandte Analysis und Stochastik, 2011.[9] T. Weiland. A discretization method for the solution of Maxwell’s equations for six-componentfields. Electronics and Communication (AEÜ), 31:116–120, 1977.[10] T. Weiland. On the unique numerical solution of Maxwellian eigenvalue problems in three dimen-sions. Particle Accelerators (PAC), 17:277–242, 1985.DOI 10.20347/WIAS.PREPRINT.2452 Berlin 2017

A multilevel Schur complement preconditioner for complex symmetric matrices

This paper describes a multilevel preconditioning technique for solving
complex symmetric sparse linear systems. The coefficient matrix is first
decoupled by domain decomposition and then an approximate inverse of the
original matrix is computed level by level. This approximate inverse is based
on low rank approximations of the local Schur complements. For this, a
symmetric singular value decomposition of a complex symmetric matix is used.
Using the example of Maxwells equations the generality of the approach is
demonstrated

Repositorium für Naturwissenschaften und Technik  (TIB Hannover)

A multilevel Schur complement preconditioner for complex symmetric
matrices

https://oa.tib.eu/renate/bitstream/123456789/2724/1/1009228048.pdf

A multilevel Schur complement preconditioner for complex symmetric matrices

Abstract

Similar works

Full text

Available Versions

Publications Server of the Weierstrass Institute for Applied Analysis and Stochastics

Repositorium für Naturwissenschaften und Technik (TIB Hannover)