Efficient cyclic reduction for QBDs with rank structured blocks: - PowerPoint PPT Presentation

Dario A. Bini, Stefano Massei, Leonardo Robol Efficient cyclic reduction for QBDs with rank structured blocks: algorithm and analysis University of Pisa Scuola Normale Superiore KU Leuven stefano.massei@sns.it Budapest, 28 June 2016 Speaker: Stefano Massei 28 June 2016 1 / 19

Quasi-birth-death processes A QBD process, in discrete time, is a bidimensional Markov chain whose transition probability matrix has the tridiagonal block Toeplitz structure   B 0 B 1 A − 1 A 0 A 1     A − 1 A 0 A 1   P = ,   ...   A − 1 A 0     ... ... with A i , B i ∈ R m × m ( m ∈ N ∪ { + ∞} ) non negative and P stochastic. m 0 0 Speaker: Stefano Massei 28 June 2016 2 / 19

The main problem Suppose m < ∞ and let the matrix P be irreducible and nonperiodic. We consider the computation of the stationary distribution of the QBD, i.e. an infinite vector π π π such that π T P = π π T , π π π π π π ≥ 0 , and � π π π � 1 = 1 . Due to the matrix-geometric property of π π , a crucial step consists in finding the π minimal non negative solution G of the quadratic matrix equation: X = A − 1 + A 0 X + A 1 X 2 , X ∈ R m × m . Many numerical methods have been proposed to address the problem and most of them are designed to deal with the general case where the block coefficients A − 1 , A 0 and A 1 have no particular structure. Speaker: Stefano Massei 28 June 2016 3 / 19

Cyclic Reduction The method on which we are going to focus is the Cyclic Reduction. Its iterative scheme requires the computation of four sequences of matrices, A ( h ) , i ( h ) , which follow the recurrence relations: i = − 1 , 0 , 1 and ˆ A 0 0 ) − 1 A ( h ) A ( h + 1 ) = A ( h ) ( I − A ( h ) 1 , 1 1 0 ) − 1 A ( h ) 0 ) − 1 A ( h ) A ( h + 1 ) = A ( h ) + A ( h ) ( I − A ( h ) − 1 + A ( h ) − 1 ( I − A ( h ) 1 , 0 0 1 0 ) − 1 A ( h ) A ( h + 1 ) = A ( h ) − 1 ( I − A ( h ) − 1 , − 1 ( h + 1 ) = ˆ ( h ) + A ( h ) 0 ) − 1 A ( h ) ˆ ( I − A ( h ) A 0 A 0 − 1 . 1 ( 0 ) = A 0 . with A ( 0 ) = A i , i = − 1 , 0 , 1 and ˆ A 0 i After each step, an approximation of the matrix G is provided by ( h ) ) − 1 A − 1 . ( I − ˆ A 0 Under mild hypothesis applicability and quadratic convergence are guaranteed. The cost of each iteration is O ( m 3 ) because it involves four matrix multiplications and the resolution of 2 m linear systems of size m . Speaker: Stefano Massei 28 June 2016 4 / 19

Cyclic Reduction/ Banded blocks For example, consider the case in which A i is finite m tridiagonal for i = − 1 , 0 , 1 (Double Quasi-Birth and Death). 0 The band structure is lost immediately when applying CR due to the inversions in its iteration scheme. Goal: Find an alternative structure to exploit for speeding up the cyclic reduction. Speaker: Stefano Massei 28 June 2016 5 / 19

Quasiseparable rank Definition A ∈ R m × m has quasiseparable rank k if the maximum rank among the off diagonal submatrices of A is k. Properties: (i) q rank ( A + B ) ≤ q rank ( A ) + q rank ( B ) (ii) q rank ( A · B ) ≤ q rank ( A ) + q rank ( B ) (iii) q rank ( A ) = q rank ( A − 1 ) Vandebril, Van Barel and Mastronardi. Matrix computations and semiseparable matrices. Johns Hopkins University Press, 2008. Speaker: Stefano Massei 28 June 2016 6 / 19

Experimental evidences Example: Cyclic reduction with starting Exponential decay of the Singular values 0 points A i ∈ R 1000 × 1000 10 −5 10 tridiagonal. Distribution of the singular −10 10 values of the sub-block −15 10 −20 10 −25 10 −30 10 −35 10 −40 10 in A ( h ) 0 5 10 15 20 25 30 35 40 during the iteration. n−th Singular value 0 Speaker: Stefano Massei 28 June 2016 7 / 19

Hierarchical matrices Strategy: Partition the row and column indices recursively and - at each step - represent the new off-diagonal blocks as low-rank outer products. Stop when the diagonal blocks reach a small dimension and represent them as full matrices. Matrix operations: Addition O ( m log ( m )) O ( m log ( m ) 2 ) Multiplication O ( m log ( m ) 2 ) Lin. System Storage O ( m log ( m )) B¨ orm, Grasedyck and Hackbusch. Hierarchical matrices. Lecture notes, 2003. Hackbusch. Hierarchical Matrices: Algorithms and Analysis. Springer Berlin, 2016. Speaker: Stefano Massei 28 June 2016 8 / 19

A bunch of applications H -matrix approximation of BEM matrices. [Hackbusch, Sauter,. . . 1990s] Matrix sign function iteration in H -arithmetic for solving matrix Lyapunov and Riccati equations. [Grasedyck,Hackbusch,Khoromskij 2004] Contour integral + H -matrices for matrix functions [Gavrilyuk et al. 2002]. H -matrix based preconditioning for FE discretization of 3D Maxwell [Ostrowski et al. 2010]. Block low-rank approximation of kernel matrices [Si, Hsieh, Dhillon 2014, Wang et al. 2015]. Clustered low-rank approximation of graphs [Savas, Dhillon 2011]. Cyclic reduction + H -matrices for quadratic matrix equations with quasiseparable coefficients. [Bini, M., Robol 2016] . . . Speaker: Stefano Massei 28 June 2016 9 / 19

Numerical Results/ Tridiagonal: Size VS Execution Time Execution time 10 5 CR H 10 − 16 H 10 − 12 10 3 H 10 − 8 Time (s) 10 1 10 − 1 10 2 10 3 10 4 10 5 10 6 Size CR H 10 − 16 H 10 − 12 H 10 − 8 Size Time (s) Residue Time (s) Residue Time (s) Residue Time (s) Residue 100 6 . 04 e − 02 1 . 91 e − 16 2 . 21 e − 01 1 . 79 e − 15 2 . 04 e − 01 8 . 26 e − 14 1 . 92 e − 01 7 . 40 e − 10 200 1 . 88 e − 01 2 . 51 e − 16 5 . 78 e − 01 1 . 39 e − 14 5 . 03 e − 01 1 . 01 e − 13 4 . 29 e − 01 2 . 29 e − 09 400 1 . 61 e + 01 2 . 09 e − 16 3 . 32 e + 00 1 . 41 e − 14 2 . 60 e + 00 1 . 33 e − 13 1 . 98 e + 00 1 . 99 e − 09 800 2 . 63 e + 01 2 . 74 e − 16 4 . 55 e + 00 1 . 94 e − 14 3 . 49 e + 00 2 . 71 e − 13 2 . 63 e + 00 2 . 69 e − 09 1600 8 . 12 e + 01 3 . 82 e − 12 1 . 18 e + 01 3 . 82 e − 12 8 . 78 e + 00 3 . 82 e − 12 6 . 24 e + 00 3 . 39 e − 09 3200 6 . 35 e + 02 5 . 46 e − 08 3 . 12 e + 01 5 . 46 e − 08 2 . 21 e + 01 5 . 46 e − 08 1 . 51 e + 01 5 . 43 e − 08 6400 5 . 03 e + 03 3 . 89 e − 08 7 . 83 e + 01 3 . 89 e − 08 5 . 38 e + 01 3 . 89 e − 08 3 . 58 e + 01 3 . 87 e − 08 12800 4 . 06 e + 04 1 . 99 e − 08 1 . 94 e + 02 1 . 99 e − 08 1 . 29 e + 02 1 . 99 e − 08 8 . 37 e + 01 1 . 97 e − 08 Speaker: Stefano Massei 28 June 2016 10 / 19

Quasiseparable rank’s growth 36 34 Q rank of A 0 32 30 Quasiseparable rank 28 26 24 22 20 18 1 2 3 4 10 10 10 10 Dimension of blocks Speaker: Stefano Massei 28 June 2016 11 / 19

Numerical Results/ Size= 1600 : Band VS Execution Time Execution time CR H 10 − 16 H 10 − 12 H 10 − 8 10 2 Time (s) 10 1 10 0 10 1 10 2 10 3 Band CR H 10 − 16 H 10 − 12 H 10 − 8 Band Time (s) Residue Time (s) Residue Time (s) Residue Time (s) Residue 2 7 . 47 e + 01 2 . 11 e − 16 1 . 58 e + 01 6 . 95 e − 15 1 . 08 e + 01 2 . 62 e − 13 7 . 86 e + 00 2 . 57 e − 09 4 7 . 65 e + 01 1 . 66 e − 16 1 . 92 e + 01 4 . 88 e − 15 1 . 48 e + 01 2 . 36 e − 13 9 . 44 e + 00 3 . 15 e − 09 8 7 . 82 e + 01 1 . 48 e − 16 2 . 81 e + 01 6 . 11 e − 15 2 . 15 e + 01 2 . 08 e − 13 1 . 31 e + 01 2 . 10 e − 09 16 7 . 50 e + 01 1 . 35 e − 16 4 . 99 e + 01 4 . 98 e − 15 3 . 48 e + 01 2 . 29 e − 13 2 . 28 e + 01 2 . 08 e − 09 32 7 . 97 e + 01 1 . 33 e − 16 9 . 40 e + 01 5 . 79 e − 15 6 . 32 e + 01 2 . 01 e − 13 4 . 15 e + 01 2 . 28 e − 09 64 8 . 03 e + 01 1 . 31 e − 16 1 . 97 e + 02 6 . 79 e − 15 1 . 29 e + 02 1 . 99 e − 13 8 . 37 e + 01 2 . 01 e − 09 128 7 . 53 e + 01 1 . 28 e − 16 4 . 01 e + 02 5 . 89 e − 15 2 . 71 e + 02 2 . 02 e − 13 1 . 75 e + 02 2 . 15 e − 09 Speaker: Stefano Massei 28 June 2016 12 / 19

Theoretical analysis/ Functional interpretation of CR We associate at each step of the CR the matrix Laurent polynomial ϕ ( h ) ( z ) := − z − 1 · A ( h ) − 1 + ( I − A ( h ) 0 ) − z · A ( h ) ϕ ( z ) := ϕ ( 0 ) ( z ) , 1 , and the matrix function defined by recurrence 2 h − 1 � ψ ( 0 ) ( z ) := ϕ ( z ) − 1 ⇒ ψ ( h ) ( z 2 h ) = 1 � ψ ( 0 ) ( ζ j z ) ψ ( h + 1 ) ( z 2 ) := 1 2 ( ψ ( h ) ( z ) + ψ ( h ) ( − z )) 2 h j = 0 Theorem (Bini,Meini) Let z ∈ C � { 0 } be such that det ( ϕ ( i ) ( z ))) � = 0 ∀ i = 0 , . . . , h and let det ( I − A ( i ) 0 ) � = 0 ∀ i = 0 , . . . , h − 1 then ϕ ( i ) ( z ) = ψ ( i ) ( z ) − 1 , i = 0 , . . . , h . � − 1 In particular ϕ ( h ) ( z 2 h ) = ψ ( h ) ( z 2 h ) − 1 = � � 2 h − 1 1 j = 0 ψ ( 0 ) ( ζ j z ) . 2 h Speaker: Stefano Massei 28 June 2016 13 / 19