A class of efficient preconditioners with multilevel sequentially semiseparable matrix structure

(1)

A Class of Efﬁcient Preconditioners with Multilevel

Sequentially Semiseparable Matrix Structure

Yue Qiu

∗

, Martin B. van Gijzen

†

_{, Jan-Willem van Wingerden}

∗

_{and Michel}

Verhaegen

∗

∗_{Delft Center for System and Control, Delft University of Technology, 2628CD Delft, the Netherlands} †_{Delft Institute of Applied Mathematics, Delft University of Technology, 2628CD Delft, the Netherlands}

Abstract. This paper presents a class of preconditioners for sparse systems arising from discretized partial differential equations (PDEs). In this class of preconditioners, we exploit the multilevel sequentially semiseparable (MSSS) structure of the system matrix. The off-diagonal blocks of MSSS matrices are of low-rank, which enables fast computations of linear complexity. In order to keep the low-rank property of the off-diagonal blocks, model reduction algorithm is necessary. We tested our preconditioners for 2D convection-diffusion equation, the computational results show the excellent performance of this approach.

Keywords: preconditioners, partial differential equations (PDEs), multilevel, sequentially semiseparable matrices PACS: 15-06, 15A06, 15A23, 15B99

INTRODUCTION

Many simulations in science and engineering design [1] [2] require the solution of a large sparse system of the form

Ax= b , (1)

where A= [ai j] is an n × n matrix and b is a given right-hand-side vector of compatible dimension. Often this is the

most time-consuming part of a simulation. If such systems originate from discretized partial differential equations (PDEs), they are usually sparse and of large size. Many efforts have been dedicated to ﬁnd efﬁcient solution methods for such systems. Some of the most popular iterative solution methods are the conjugate gradient (CG), minimal residual (MINRES), generalized minimal residual (GMRES) and induced dimension reduction (IDR(s)) methods [3] [4] [5]. The performance of these methods highly depends on the choice of the preconditioners. In this paper, we study a new class of preconditioners based on the multilevel sequentially semiseparable matrix structure of the system [6].

Semiseparable matrices appear in several types of applications, e.g. the ﬁeld of integral equations [7], Gauss-Markov processes [8], boundary value problems [9] and rational interpolation [10]. Semiseparable matrices are deﬁned in [11] as matrices of which all sub-matrices taken from the lower-triangular or upper-triangular part are of rank 1 . The generalization of semiseparable matrices, called quasiseparable matrices, are matrices of which all sub-matrices extracted from the strictly lower-triangular or strictly upper-triangular part are of rank 1 [11], where semiseparable plus diagonal matrices are of such type. This more general type of matrices is called sequentially semi-separable (SSS) by Dewilde et al [12]. For this class, the off-diagonal blocks are of low-rank, but not necessarily limited to 1. The property of low-rank off-diagonal blocks is also investigated by Eidelman et al [13]. They also call this general class of matrices quasiseparable matrices. The class of SSS matrices is closed under basic operations as addition, multiplication and inversion. Moreover, decompositions/factorizations such as the QR [14], LU/LDU [11] [15], and ULV factorizations [16] [17] can also be computed in a structure preserving way. Many operations on SSS matrices can be performed with linear computational complexity. Examples are the matrix-matrix multiplication, the matrix-vector multiplication, matrix inversion, and the QR, LU and ULV factorizations [18] [16] [14].

Systems that arise from the discretisation of 1D partial differential equations typically have an SSS structure. Discretisations of higher dimensional (2D or 3D) partial differential equations give rise to matrices that have a so-called multi-level SSS structure [15]. For such MSSS matrices, operations like inversion and decomposition can not be performed with linear complexity. However, as discussed in [15], and as we will show in this paper, it is possible to compute an inexact LU decomposition of an MSSS matrix that can be used as a preconditioner. By numerical experiments we will show that this preconditioner can be computed at linear complexity and that the number of iterations is almost independent of the mesh size.

11th International Conference of Numerical Analysis and Applied Mathematics 2013

2253

This article is copyrighted as indicated in the abstract. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 131.180.131.94 On: Fri, 25 Oct 2013 09:43:41

(2)

The outline of this paper is as follows. Section 2 introduces multilevel sequentially semiseparable matrices and explains their relation with discretized PDEs. Section 3 discussed the MSSS preconditioner and illustrates its perfor-mance on two convection-diffusion problems. The last section gives the conclusions.

MULTILEVEL SEQUENTIALLY SEMISEPARABLE MATRICES

The matrices in this paper will always be real and their dimensions are compatible for matrix- matrix addition, multiplication and matrix-vector multiplication when their sizes are not mentioned. To keep this paper self-contained, we review some definitions and concepts, see also [6]. The generators representation of sequentially semiseparable matrices are defined by Definition 1.

Deﬁnition 1. Let A be an N × N matrix with the SSS structure. Let m1, m2, ··· mn be positive integers with

N= m1+ m2+ ··· + mnsuch that A can be written in the following block-partitioned form:

Ai j= ⎧ ⎨ ⎩ UiWi+1···Wj−1VjT, i < j; Di, i= j; PiRi−1···Rj+1QTj, i> j (2) where superscript ’T ’ denotes the transpose of the matrix. The above representation of A is called the generators representation. The sequences{Ui}n−1i=1,{Wi}in−1=2,{Vi}ni=2,{Di}ni=1,{Pi}ni=2,{Ri}n−1i=2,{Qi}n−1i=1 are matrices whose

sizes are listed in table 1 and they are called generators of the SSS matrix A. TABLE 1. Size of the generators of A

Generators Ui Wi Vi Di Pi Ri Qi

Sizes mi× ki ki−1× ki mi× ki−1 mi× mi mi× li li−1× li mi× li+1

With the generators parametrization, the SSS matrix A can be expressed by its generators and denoted as

A= S S S (Ps,Rs,Qs,Ds,Us,Ws,Vs). (3)

Take n= 4 for example, the SSS matrix A is shown by (4), ⎡ ⎢ ⎢ ⎣ D1 U1V2T U1W2V3T U1W2W3V4T P2QT1 D2 U2V3T U2W3V4T P3R2QT1 P3QT2 D3 U3V4T P4R3R2QT1 P4R3QT2 P4QT3 D4 ⎤ ⎥ ⎥ ⎦ (4)

Using this representation, the structure of SSS matrices can be exploited to do fast computation for matrix addition, multiplication and inversion by operations on its generators, with the cost of linear computational complexity. Table 2 lists references to papers that discuss how a given matrix operation can be performed using SSS matrix arithmetic.

TABLE 2. Commonly used operations on SSS matrices

operations A∗ x A± B A∗ B A−1 LU Model Reduction Lx= b∗

references [12] [13] [18] [12] [13] [18] [12] [13] [18] [13] [6] [11] [15] [6] [18] [18]

∗_{L is a lower triangular SSS matrix}

Remark 1. With the operations listed in table 2, it is always possible to construct a preconditioner and solve the

preconditioned system with linear complexity.

Analogously to SSS matrices are deﬁned in [15] by Deﬁnition 2.

Deﬁnition 2. The matrix A is said to be a k-level MSSS matrix where all its generators are (k −1)-level MSSS matrices.

The 1-level MSSS matrix is the SSS matrix that satisﬁes Deﬁnition 1.

2254

(3)

Example 1. Any banded matrix is automatically semiseparable according to [11]. Similarly, every block banded

matrix is multilevel sequentially semiseparable (MSSS). Take for example the discretized 2D Laplace equation with homogeneous Dirichlet boundary conditions by ﬁnite-element discretization on a regular grid, the stiffness matrix K

is given by K= ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ D F F D F F . .. ... . .. F F D ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , where D= ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 4 −1 −1 4 −1 −1 . .. . .. −1 −1 4 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦

and F= −I where I is the

identity matrix. The matrices D and F have the SSS matrix structure and can be denoted as

D= S S S (1, 0, −1, 4, 1, 0, −1) and F = S S S (0, 0, 0, −1, 0, 0, 0). (5) Meanwhile, the matrix K has the (2-level) MSSS matrix structure and can be denoted by

K= M S S S (I, 0, F, D, I, 0, F). (6)

MSSS PRECONDITIONERS

In this section we will describe how the MSSS structure of a matrix can be exploited to construct an efﬁcient preconditioner. We will illustrate the approach with experiments on a convection-diffusion problem.

To construct a preconditioner, we factorize the system matrix with MSSS matrix techniques. The decompo-sition method proposed in [15] results in a growth of the rank of the generators. As a result, the decomposi-tion of an MSSS matrix is not of order N complexity. The soludecomposi-tion proposed in [15] is to use a model or-der reduction algorithm to reduce the rank of the sub-blocks, to keep the decomposition operation of oror-der N. However, as a results of the model order reduction, the decomposition is only approximate. The model reduc-tion algorithm plays a key rule in MSSS arithmetic. In [15] Gondzio et al used the method in [19] to solve the PDE-constrained optimization problem. We use a novel model reduction algorithm described in [6] (available at https://www.dropbox.com/sh/pj6gqw6pxzn69qz/TXyboPQEhy) to construct an efficient preconditioner for the 2D convection-diffusion problem. The idea of the novel model reduction algorithm is to use the low-rank property of the controllability and observability gramians of linear time-varying (LTV) systems to approximate the gramians with low numerical rank. With the approximate gramians of LTV systems, it is possible to perform the balanced truncation to get a LTV system of reduced order, which corresponds to the reduced generators of an SSS matrix. The model reduction in [19] is to approximate the Hankel blocks of an SSS matrix. Compared with this model reduction algo-rithm, the novel model reduction algorithm needs less floating point operations. Meanwhile, the novel model reduction algorithm leads to satisfied accuracy. Numerical experiments show the efficiency of our method.

Example 2. Let Ω = [0, 1]2and consider the problem

−ε∇2_u_{+ −}→_{ω .∇u = 0 in Ω} ₍₇₎ with u= uDonΩDand _∂n∂u= uNonΩNwhere uD= uN= (2x − 1)2(2y − 1)2,ΩD= {0} × [0, 12], ΩN= [0, 12] × {0} andε is a positive scalar, −→ω is the unit directional vector, −→ω = (cos(θ), sin(θ))T.

The discretized 2D convection-diffusion equation yields a nonsymmetric linear system, which we solve with the induced dimension reduction (IDR(s)) method in combination with the MSSS preconditoner. The experiments have been performed on a laptop with Intel Core 2 Duo 2.4 GHZ and 4Gb memory using Matlab R2010b. The parameters are set toθ =π₅, s= 4, the maximum semiseparable order is listed in the brackets followed after the mesh size and the computational results for differentε are shown in table 3 and 4

Remark 2. Table 3 and 4 show that the number of iterations is almost independent of the mesh size, and that the

computing time of both the computation of the preconditioner and of the iterative solution increases linearly with the number of equations.

2255

(4)

TABLE 3. Number of iterations and time of preconditioning and IDR(s) solver whenε = 10−1 mesh size h iterations Preconditioning (seconds) IDR(s) (seconds) total (seconds)

2−5(3) 4 0.56 0.20 0.76

2−6(3) 4 1.91 0.36 2.27

2−7(4) 4 6.34 0.77 7.11

2−8(4) 4 23.91 2.09 26.00

2−9(4) 5 92.08 9.95 102.03

TABLE 4. Number of iterations and time of preconditioning and IDR(s) solver whenε = 10−2 mesh size h iterations Preconditioning (seconds) IDR(s) (seconds) total (seconds)

2−5(3) 2 0.65 0.16 0.81 2−6(3) 3 1.99 0.29 2.28 2−7(3) 3 5.97 0.67 6.64 2−8(4) 4 23.27 1.71 24.98 2−9(4) 4 90.05 6.36 96.41

CONCLUSIONS

In this paper, we have exploited the MSSS structure of the system matrix to construct an efficient preconditioner. The MSSS structure has low-rank off-diagonal blocks, which enables fast computations of linear complexity. In order to keep the off-diagonal blocks of low-rank, a novel model reduction algorithm is used. Numerical experiment of 2D convection-diffusion problem shows the efficiency of this approach. This approach can be extended to preconditioning for more complicated 3D problems. Moreover, our method is also efficient for preconditioning of PDE-constrained optimization problems, we refer [6] for further results.

REFERENCES

1. L. Biegler, O. Ghattas, M. Heinkenschloss, D. Keyes, and B. van Bloemen Waanders, Real-time PDE-constrained Optimization, Society for Industrial and Applied Mathematics, 2007.

2. T. Rees, H. Dollar, and A. Wathen, SIAM Journal on Scientiﬁc Computing 32, 271–298 (2010). 3. G. Golub, and C. Van Loan, Matrix computations, Johns Hopkins University Press, 1996.

4. Y. Saad, Iterative methods for sparse linear systems, Society for Industrial and Applied Mathematics, 2003. 5. P. Sonneveld, and M. van Gijzen, SIAM Journal on Scientiﬁc Computing 31, 1035–1062 (2008).

6. Y. Qiu, M. van Gijzen, J. van Wingerden, and M. Verhaegen, Efﬁcient preconditioners for PDE-constrained optimization problems with a multi-level sequentially semi-separable matrix structure, Tech. Rep. 13-04, Delft Institution of Applied Mathematics, Delft University of Technology (2013), URL https://www.dropbox.com/sh/pj6gqw6pxzn69qz/ TXyboPQEhy.

7. R. Gonzales, J. Eisert, I. Koltracht, M. Neumann, and G. Rawitscher, Journal of Computational Physics 134, 134–149 (1997). 8. A. Kavcic, and J. Moura, IEEE Transactions on Information Theory 46, 1495–1509 (2000), ISSN 0018-9448.

9. L. Greengard, and V. Rokhlin, Communications on Pure and Applied Mathematics 44, 419–452 (1991).

10. M. Van Barel, D. Fasino, L. Gemignani, and N. Mastronardi, “Orthogonal rational functions and diagonal-plus-semiseparable matrices,” in International Symposium on Optical Science and Technology, 2002, pp. 162–170.

11. R. Vandebril, M. Van Barel, and N. Mastronardi, Matrix computations and semiseparable matrices: linear systems, Johns Hopkins University Press, 2007.

12. S. Chandrasekaran, P. Dewilde, M. Gu, T. Pals, X. Sun, A. van der Veen, and D. White, SIAM Journal on Matrix Analysis and Applications27, 341–364 (2005).

13. Y. Eidelman, and I. Gohberg, Integral Equations and Operator Theory 34, 293–324 (1999).

14. Y. Eidelman, I. Gohberg, and V. Olshevsky, Linear Algebra and Its Applications 404, 305–324 (2005). 15. J. Gondzio, and P. Zhlobich, arXiv preprint arXiv:1112.6018 (2011).

16. S. Chandrasekaran, M. Gu, and T. Pals, SIAM Journal on Matrix Analysis and Applications 28, 603–622 (2006).

17. A. van der Veen, Time-Varying System Theory and Computational Modeling, Ph.D. thesis, Delft University of Technology (1993).

18. S. Chandrasekaran, P. Dewilde, M. Gu, T. Pals, and A. van der Veen, Fast stable solvers for sequentially semi-separable linear systems of equations, Tech. rep., Lawrence Livermore National Laboratory (2003).

19. Y. Eidelman, and I. Gohberg, Calcolo 42, 187–214 (2005).

2256