Efficient cyclic reduction for QBDs with rank structured blocks∗ Dario A. Bini (bini@dm.unipi.it), Dipartimento di Matematica, University of Pisa, Italy Stefano Massei (stefano.massei@sns.it), Scuola Normale Superiore, Pisa, Italy Leonardo Robol (leonardo.robol@sns.it), 6 1 Scuola Normale Superiore, Pisa, Italy 0 2 January 6, 2016 n a J 5 Abstract ] WeprovideeffectivealgorithmsforsolvingblocktridiagonalblockToeplitzsystemswith A m×mquasiseparableblocks,aswellasquadraticmatrixequationswithm×mquasiseparable N coefficients, based on cyclic reduction and on the technology of rank-structured matrices. The algorithms rely on the exponential decay of the singular values of the off-diagonal . h submatrices generated by cyclic reduction. We provide a formal proof of this decay in the t a Markovian framework. The results of the numerical experiments that we report confirm a m significantspeedupoverthegeneralalgorithms,alreadystartingwiththemoderatelysmall size m≈102. [ 1 v 1 Introduction 1 6 Cyclic reduction (CR) is an effective tool that can be used for solving several problems in linear 8 algebra and in polynomial computations [4]. It has been originally introduced by G.H. Golub 0 andR.W.Hockneyinthemid1960s[13,8],forthenumericalsolutionofblocktridiagonallinear 0 . systems stemming from the finite differences solution of elliptic problems, and has been gener- 1 alized to solve nonlinear matrix equations associated with matrix power series with applications 0 6 to queuing problems, Markov chains and spectral decomposition of polynomials. In fact, an 1 important application of CR concerns the computation of the minimal nonnegative solution of : thematrixequationX =A +A X+A X2,encounteredinQuasiBirth-Death(QBD)Markov v −1 0 1 i chains, where A−1, A0, and A1 are given m×m nonnegative matrices such that A−1+A0+A1 X is irreducible and stochastic and where X is the matrix unknown [2, 4]. The computation of the r solution X allows to recover the steady state vector π of the Markov chain. a RewritingthematrixequationasA +(A −I)X+A X2 =0,CRcomputesfoursequences −1 0 1 ∗Work supported by Gruppo Nazionale di Calcolo Scientifico (GNCS) of INdAM, and by the PRA project “Mathematicalmodelsandcomputationalmethodsforcomplexnetworks”fundedbytheUniversityofPisa. 1 of matrices, A(h), i=−1,0,1 and A(cid:98)(h), according to the following equations i 0 A(h+1) =−A(h) S(h) A(h), S(h) =(A(h)−I)−1 1 1 1 0 A(h+1) =A(h)−A(h) S(h) A(h) − A(h) S(h) A(h), (1) 0 0 1 −1 −1 1 A(h+1) =−A(h) S(h) A(h), A(cid:98)(h+1) =A(cid:98)(h)−A(h) S(h) A(h), −1 −1 −1 0 0 1 −1 for h=0,1,..., with A(i0) =Ai, i=−1,0,1 and A(cid:98)(00) =A0−I. It can be proved, in the context of Markov chains, [2] that the sequence −(A(cid:98)0(h))−1A−1 converges to the minimal nonnegative solution G of the matrix equation. If the QBD process is not null recurrent, applicability and quadratic convergence of the algorithm are guaranteed [2][Theorems 7.5, 7.6]. In the null recurrent case, it has been proved in [11] that convergence is linear with factor 1. 2 Without any further assumption on the structure of the blocks, each step of CR requires a smallnumberofmatrixmultiplicationsandonematrixinversionfortheresultingcomputational costof O(m3)arithmeticoperations(ops)perstep. Ontheotherhand, thereareseveralmodels from the applications in which the blocks A exhibit special structures. In order to decrease the i computational complexity of the iterations, variations of CR which exploit these structures have been proposed, see for instance [3], [17], [1]. Here,weareinterestedinanalysingthecasewheretheblocksA arequasiseparablematrices. i That is, the case where the off-diagonal submatrices of A , A and A , strictly contained in −1 0 1 the upper or in the lower triangular part, have low rank with respect to m. The maximum of the ranks of the off-diagonal submatrices is called quasiseparable rank and a matrix with quasiseparablerankk iscalledk-quasiseparable. Observethatk-quasiseparablematricesinclude bandedmatrices. Thesestructuresareencounteredinwideandimportantclassesofapplications like, bidimensional random walks [16], the Jackson tandem queue model [14] and other QBD processes, or, for instance, in the finite differences discretization of elliptic PDEs. Our goal is to design a version of CR which exploits the rank structures of the blocks A and i which can be implemented at a substantially lower cost. This way, we may arrive at designing effective solvers both for block tridiagonal block Toeplitz systems and for the quadratic matrix equations encountered in QBD Markov chains. Indeed, a way to reach this goal is to find out if some structure of the blocks A(h) is maintained during the CR steps (1), and then to try to i exploit this structure in order to design an effective implementation of CR. Looking at the iterative scheme (1), one can find out that the quasiseparable rank can grow exponentially at each step. Despite that, plotting the singular values of the off-diagonal blocks of the matrices A(h) shows an interesting behaviour as reported in Figure 1. For a randomly i generatedQBDprocesswithtridiagonalblocksofsize1600,thesingularvaluesofanoff-diagonal submatrix of size 799×800 in A(h) have an exponential decay in the first 20 steps of CR. 0 It is evident that, even though the number of nonzero singular values grows at each step of CR, the number of singular values above the machine precision – denoted by a horizontal line in Figure 1– is bounded by a moderate constant. Moreover, the singular values seem to stay below a straight-line which constitutes an asymptotic bound. That is, they get closer to this line as h→∞. Thelogarithmscalesuggestthatthecomputedsingularvaluesσ(h) decayexponentially i with i and the basis of the exponential grows with h but has a limit less than 1. Wewillprovethisasymptoticpropertyrelyingonthetechnologyofrank-structuredmatrices and relate the basis of the exponential decay to the width of the domain of analiticity of the matrix function ϕ(z)−1 for ϕ(z)=−z−1A +I−A −zA . −1 0 1 The paper is organized as follows. In Section 2 we recall some useful properties of CR. In Section 3 we prove some preliminary lemmas for bounding the singular values of sums, products 2 100 Exponential decay of the Singular values 10−5 10−10 10−15 10−20 10−25 10−30 10−35 10−400 5 10 15 20 25 30 35 40 n−th Singular value Figure 1: Log-scale plot of the largest singular values of the largest south-western submatrix of A containedinthelowertriangularpart,form=1600. Thehorizontallinedenotesthemachine 0 precision threshold. Matrices are randomly generated so that A (cid:62) 0 are tridiagonal matrices i and A +A +I+A is stochastic. −1 0 1 andinversesofmatrices. Section4containstheproofofthemainresultconcerningtheexponen- tial decay of the singular values. In Section 5 we describe our algorithm which implements CR relying on the package H2Lib [5], concerning H-matrices, and the related code together with the results of the numerical experiments. Finally, in Section 6, we give some concluding comments and possible future developments. 2 Main properties of cyclic reduction We recall a functional interpretation of cyclic reduction introduced in the Markov chains framework [2], [4] to prove applicability and convergence properties. Associate the matrices A(h), i = −1,0,1 defined in (1) with the matrix Laurent polynomial i ϕ(h)(z) := −z−1A(h) +(I −A(h))−zA(h), where ϕ(0)(z) = ϕ(z) = z−1A +(I −A )−zA , −1 0 1 −1 0 1 and define the matrix rational function ψ(h)(z)=ϕ(h)(z)−1. The following property holds (cid:40) ψ(0)(z):=ψ(z), ψ(h+1)(z2):= 1(ψ(h)(z)+ψ(h)(−z)). 2 In particular, making explicit the recurrence relation in the sequence {ψ(h)}, we find that ψ(h)(z2h) = 21h (cid:80)2j=h−01ψ(0)(ωNj z) where ωN = e2Nπi is a principal N-th root of the unity for N =2h, and i denotes the imaginary unit, so that  −1 2h−1 ϕ(h)(z2h)= 1 (cid:88) ψ(0)(ωj z) . (2) 2h N j=0 Observe that in the case where A , A and A are tridiagonal, then ϕ(z) is tridiagonal as −1 0 1 well, so that for any value of z such that detϕ(z)(cid:54)=0, the matrix ψ(z) is semiseparable, that is, tril(ψ(z))=tril(L), triu(ψ(z))=triu(U), where L and U are matrices of rank 1 [18]. 3 Relying on this property, in the next section we will show that ψ (z), as well as the blocks h A(h),A(h) andA(h),hasoff-diagonalsubmatriceswithsingularvalueswhichdecayexponentially. −1 0 1 3 Singular values of sums, products and inverses Inthissection,werecallsomebasicfactsonthesingularvaluesdecomposition(SVD)andprovide some technical lemmas. Let us denote by σ (A) the j-th singular value of the m×n matrix A j sorted in non-ascending order and by A = UΣV∗, Σ = diag(σ ,...,σ ), p = min{m,n}, the 1 p SVD of A where U and V are unitary and V∗ denotes the transpose conjugate of V. Moreover, denotebyu andv theithcolumnofU andV, respectively. Werecallthefollowingwell-known i i property concerning the SVD of a matrix A. Property 3.1. For any matrix A the function f(X) = (cid:107)A−X(cid:107) takes its minimum over the 2 class of matrices of rank (cid:96)−1 for X =(cid:80)(cid:96)−1σ u v∗ and the value of the minimum is exactly σ . i=1 i i i (cid:96) 3.1 Some technical lemmas Let (cid:107)A(cid:107) be the Euclidean norm of A and µ (A) = (cid:107)A(cid:107) ·(cid:107)A−1(cid:107) be the condition number of 2 2 2 2 A. The following result relates the singular values of the matrices A, B and C =AB. Lemma 3.2. Consider two matrices A∈Cm×n, B ∈Cn×n, such that B is invertible. Then σ (A) σ (A) j (cid:54)σ (AB)(cid:54)(cid:107)B(cid:107) ·σ (A), j (cid:54)σ (BA∗)(cid:54)(cid:107)B(cid:107) ·σ (A), (cid:107)B−1(cid:107) j 2 j (cid:107)B−1(cid:107) j 2 j 2 2 Proof. We prove only the first statement since the other follows directly from it. Consider the SVD decompositions A = U Σ V∗, B = U Σ V∗. Recall that the singular values of a generic A A A B B B matrix M ∈ Cm×n are the square roots of the eigenvalues of M∗M, so that by the minimax theorem, we can write x∗M∗Mx σ (M)2 = max min . j dim(V)=j x∈V x∗x V⊆Cn x(cid:54)=0 Now note that AB = U Σ V∗U Σ V∗ and since unitary matrices do not change the singular A A A B B values we have σ (AB)=σ (Σ QΣ ) where Q=V∗U . Thus, we can express σ (AB)2 as j j A B A B j x∗Σ∗Q∗Σ∗Σ QΣ x (Σ x)∗Q∗Σ2Q(Σ x) σ (AB)2 = max min B A A B = max min B A B . j dim(V)=jx∈V x∗x dim(V)=jx∈V x∗x V⊆Cn V⊆Cn By setting y =Σ x and recalling that Σ is invertible by hypothesis we have B B y∗Q∗Σ2Qy y∗y σ (AB)2 = max min A · j dim(V)=jy∈V y∗y x∗x V⊆Cn so that, by using the fact that Q is unitary and that x∗x (cid:54) y∗y (cid:54) (cid:107)B(cid:107)2x∗x, we obtain (cid:107)B−1(cid:107)2 2 2 σj(A)2 = σj(ΣA)2 (cid:54)σ (AB)2 (cid:54)(cid:107)B(cid:107)2·σ (Σ )2 =(cid:107)B(cid:107)2·σ (A)2. (cid:107)B−1(cid:107)2 (cid:107)B−1(cid:107)2 j 2 j A 2 j 2 2 The following lemma relates the singular values of an off-diagonal submatrix of the inverse of a given matrix A with those of the corresponding submatrix of A. 4 Lemma 3.3. Let A∈Cn×n be an invertible matrix and consider the block partitioning (cid:18) A B (cid:19) (cid:18) A˜ B˜ (cid:19) A= , A−1 = , C D C˜ D˜ where A and A˜ are i×i matrices for 2(cid:54)i(cid:54)n−1. We have the following properties 1. If D is invertible then 1 σ (C) (cid:54) σ (C˜) (cid:54) (cid:107)D−1(cid:107) ·(cid:107)S−1(cid:107) ·σ (C), where S = (cid:107)D(cid:107)2(cid:107)SD(cid:107)2 j j 2 D 2 j D A−BD−1C is the Schur complement of D. 2. If A is invertible then 1 σ (C) (cid:54) σ (C˜) (cid:54) (cid:107)A−1(cid:107) ·(cid:107)S−1(cid:107) ·σ (C), where S = (cid:107)A(cid:107)2(cid:107)SA(cid:107)2 j j 2 A 2 j A D−CA−1B is the Schur complement of A. Proof. Let us consider part 1. If D is invertible we can write (cid:18) A B (cid:19)−1 (cid:18) S−1 S−1BD−1 (cid:19) = D D C D −D−1CS−1 D−1+D−1CS−1BD−1 D D and in particular we have C˜ = −D−1CS−1. Repeatedly applying Lemma 3.2 to C˜ yields D 1 σ (CS−1) (cid:54) σ (C˜) (cid:54) (cid:107)D−1(cid:107) · σ (CS−1) and eventually 1 σ (C) (cid:54) σ (CS−1) (cid:54) (cid:107)D(cid:107)2 j D j 2 j D (cid:107)SD(cid:107)2 j j D (cid:107)S−1(cid:107) ·σ (C). Combining these inequalities gives us the thesis. For proving part 2, we can D 2 j proceed in the same manner relying on the inversion formula (cid:18) A B (cid:19)−1 (cid:18) A−1+A−1BS−1CA−1 −A−1BS−1 (cid:19) = A A . C D −S−1CA−1 S−1 A A Lemma 3.4. Let A = (cid:80)+∞ A and A+ = (cid:80)+∞A where A ∈ Cm×n have rank at most k i=−∞ i i=0 i i and suppose that (cid:107)Ai(cid:107)2 (cid:54)Me−α|i|. Then σj(A)(cid:54) 1−2Me−α ·e−αj2−kk, σj(A+)(cid:54) 1−Me−α ·e−αj−kk. Proof. Note that B = (cid:80)s−1 A has rank at most k(2s − 1). For j positive integer, set s i=−s+1 i s=(cid:100)j−k(cid:101) and observe that, since 2k(cid:100)j−k(cid:101) 0, let ω = ei2π/N and observe that the quantity 1 (cid:80)N−1(zωj )k N N j=0 N coincides with zk if k ≡ 0 mod N and with 0 otherwise. This way, if A(z) = (cid:80) ziA i∈Z i is a matrix Laurent series analytic on the annulus A(r ,r ) = {z ∈ C : r < |z| < r } for 1 2 1 2 5 0 < r < 1 < r , then 1 (cid:80)N−1A(ωjz) = B(zN) where B(z) = (cid:80) ziA is analytic on 1 2 N j=0 n i∈Z Ni A(rN,rN). WedenotebyI theoperatorwhichmapsA(z)intoB(z)andwriteB(z)=I (A(z)). 1 2 N N Observe that I is linear and continuous on the space of analytic functions on A(r ,r ). N 1 2 Moreover,inviewof (2),wehaveψ(h) =I (ψ(0))forN =2h. Thisway,ifweprovethatany N off-diagonalsubmatrixB(z)ofψ(0)(z)issuchthatI (B(z))hastheexponentialdecayproperty N for its singular values, then we have shown this property also for ψ(h)(z). Partition ϕ(z) and ϕ(z)−1 as follows (cid:20)I−E(z) −B(z) (cid:21) (cid:20)E˜(z) B˜(z)(cid:21) ϕ(z)= , ψ(z):=ϕ(z)−1 = , (3) −C(z) I−D(z) C˜(z) D˜(z) where E(z) and D(z) are square matrices of any compatible size. Theorem 4.1. Let ϕ(z)=−z−1A +I−A −zA be an m×m matrix function such that −1 0 1 (i) The A are non-negative and I−ϕ(z) has spectral radius smaller than 1 for any z ∈T. i (ii) The blocks A are k-quasiseparable and (cid:107)A (cid:107) (cid:54)L, i=−1,0,1. i i 2 (iii) There exist t>1 and δ >0 such that detϕ(z)(cid:54)=0 and (cid:107)ϕ(z)−1(cid:107) (cid:54)δ for z ∈A(t−1,t). 2 Then ρ(I − ϕ(z)) < 1 for any z ∈ A and in the partitioning (3), both blocks I − E(z) and I −D(z) are invertible for any z ∈ A. Moreover, for any z ∈ T and for any h, the singular values of C˜(h) :=I (C˜(z)), with N =2h, are such that N 4Lδ2 σs(C˜(h)(z))(cid:54)3Me−s−6k3klogt, M = (1−e−Nlogt)(1−t−1). (4) Moreover, if A , A , A are tridiagonal then the above bound turns into −1 0 1 σs(C˜(h)(z))(cid:54)Me−2slogt. (5) Proof. Let us prove that ρ(I −ϕ(z)) < 1 for any z ∈ A. By contradiction, assume that there exists ξ ∈A such that ρ(I−ϕ(ξ))(cid:62)1. Since A (cid:62)0 for i=−1,0,1 then |I−ϕ(ξ)|(cid:54)I−ϕ(|ξ|), i and by the monotonicity of the spectral radius we get 1 (cid:54) ρ(I −ϕ(ξ)) (cid:54) ρ(I −ϕ(|ξ|)). Thus, sinceρ(I−ϕ(1))<1(cid:54)ρ(I−ϕ(|ξ|))andρisacontinuousfunction,thenthereexists1/t<ξˆn. This allows us to write j i k ∞ (cid:32) k (cid:33) ∞ A= 1 (cid:88)A =(cid:88) 1 (cid:88)σ (A )u v∗ =(cid:88)A˜ . k i k j i i,j i,j j i=1 j=1 i=1 j=1 Observe that A˜ have rank k and (cid:107)A (cid:107)(cid:54)γe−αj. Applying Lemma 3.4 completes the proof. j j We may conclude with the decay property for the singular values of the off-diagonal submatrices of A(h), for i=−1,0,1. i Lemma 4.7. Let ϕ(h)(z) be the matrix function genertaed at the hth step of CR with the property that every offdiagonal submatrix B(z) of ϕ(h)(z) has decaying singular values such that σ (B(z)) (cid:54) γe−αs. Then every coefficient B of B(z) = z−1B + B + zB is such that s i −1 0 1 σs(B0)(cid:54)γe−αj−22, σs(Bi)(cid:54)γe−αj−33, for i=1,−1. Proof. By Lemma 4.5, we have an expression for B based on evaluations of B(z). In particular, i wehaveA = 1(B(i)+B(−i)), A = 1(ξ∓1B(ξ)+ξ∓5B(ξ5)−B(−1)),whereξ isaprimitive 0 2 ±1 3 6-th root of the unity. Applying Lemma 4.6 completes the proof. 4.3 The Markovian case One of the most interesting application of CR algorithm is in the Markovian framework, in whichapplicabilityandconvergencepropertiesareguaranteed. Inthatcase,thematrixfunction ϕ satisfies almost all the hypotheses made in the previous subsections but it is singular at z =1 since 1 is always an eigenvalue of ϕ(z). Nevertheless we will show that Corollary 4.4 can still be applied considering a rescaled version of ϕ(z). 9 When the coefficients A for i = −1,0,1 represent the blocks of the transition matrix of an i irreduciblenotnullrecurrentQBDprocess,theeigenvaluesofϕ(z)enjoythefollowingproperties [9, 4]: (i) |λ | (cid:54) |λ | (cid:54) ... (cid:54) |λ | (cid:54) λ < λ (cid:54) |λ | (cid:54) ... (cid:54) |λ |, with λ ,λ ∈ R and 1 2 m−1 m m+1 m+2 2m m m+1 one of the two equal to 1. (ii) In the annulus {λ < |z| < λ } ϕ is invertible and the spectral radius of I −ϕ(z) is m m+1 strictly less than 1. Hence we consider the rescaled version of ϕ, that is, ϕ (z) := ϕ(θz), and we choose θ = θ (cid:113) (cid:112)λ λ . We obtain a matrix function invertible on A={1 <|z| 1 and δ (cid:62) 0, consider the following class of matrix functions associated with QBD stochastic processes with k-quasiseparable blocks: χ :=(cid:8)ϕ(z): (cid:107)ϕ−1(z)(cid:107) (cid:54)δ t−1 (cid:54)|z|(cid:54)t, t<λ /λ (cid:9). δ,t 2 m+1 m Thenthereexistsauniformconstantγ(δ,t)suchthatforanyoff-diagonalblockC(h)(z)ofϕ(h)(z), with ϕ∈χδ,t, its s-th singular value is bounded by σs(C(h)(z))(cid:54)γ(δ,t)·e−s−6k3klogt. Remark4.9. Observethatinthecaseofnull-recurrentQBDprocessesonehasλ =λ =1, m m+1 so that there is no open annulus including T where ϕ(z) is nonsingular and we cannot apply Theorem 4.1. This drawback can be partially overcome by applying the shift technique of [2, 4]. This technique allows to construct a new matrix function ϕ˜(z) which has the same eigenvalues of ϕ(z) except for the eigenvalue 1 which is shifted to 0. So that ϕ˜(z) has an open annulus containingTwhereitisnonsingular. Moreover,applyingCRtoϕ˜(z)generatesmatrixsequences which easily allow to recover the corresponding matrix sequences obtained by applying CR to ϕ(z). The sequences associated with ϕ(z) differ from the sequences associated with ϕ˜(z) by a rank-1 correction. This way, if the exponential decay of the singular values holds for the latter sequences, it holds also for the former ones. The difficulty that still remains is that the nonnegativityoftheblocksA ,A andA isnotgenerallysatisfiedbythefunctionϕ˜(z)sothat −1 0 1 in principle Theorem 4.1 cannot be applied and a different version specific for this case should be formulated. Remark 4.10. The bounds that we have given to the decay of the singular values of the off- diagonalsubmatricesarenotstrict. Experimentally,singularvaluesseemtodecayslightlyfaster. Even in the null-recurrent case where λ = λ = 1, the decay still occurs even though in a m m+1 deterioratedform. Thedecaypropertiesclearlydependonthedomainofanalyticityofψ(z)but this is not the only reason of the decay. More investigation is needed in this direction. 5 An algorithm using H-matrices We have provided an implementation of CR, which applies to matrix functions ϕ(z) having quasiseparable blocks, and relies on the approximate quasiseparable structure induced by the decay of the singular values. We relied on the H-matrix representation of [6, 7, 10]. 10

Efficient cyclic reduction for QBDs with rank structured blocks

Similar Efficient cyclic reduction for QBDs with rank structured blocks