A FAST METHOD FOR SOLVING A BLOCK TRIDIAGONAL QUASI-TOEPLITZ LINEAR SYSTEM

. This paper addresses the problem of solving block tridiagonal quasi-Toeplitz linear systems. Inspired by [9], we propose a more generalized algorithm for such systems. The algorithm is based on a block decomposition for a block tridiagonal quasi-Toeplitz matrix and the Sherman-Morrison-Woodbury inversion formula. We also compare the proposed approach to the standard block LU decomposition method. A theoretical accuracy and error analysis is also considered. All algorithms have been implemented in Matlab. Numerical experiments performed with a wide variety of test problems show the eﬀectiveness of our algorithm in terms of eﬃcience, stability and robustness.


1.
Introduction. An n × n matrix A is said to be quasi Toeplitz if it is a small rank perturbation of a Toeplitz matrix.
Quasi-Toeplitz matrices arise in many mathematical and engineering investigations. For instance, the inverse of a Toeplitz matrix is quasi Toeplitz [13]. And, when we use cubic B-splines collocation method for solving nonlinear parabolic differential equations, we find that the collocation matrices are quasi-Toeplitz. Especially, when considering the Neumanns or Dirichlets type boundary conditions, the collocation matrices are tridiagonal quasi-Toeplitz [21,16,22,20]. Also, in curve or surface fitting, which are encountered frequently in CAD/CAM, the collocation matrices resulted by cubic B-splines basis are also tridiagonal quasi-Toeplitz type [15]. Block quasi Toeplitz matrices arise in the numerical approximation of time-dependent partial differential equations (PDEs) by generalizations of implicit multistep formulas used in boundary value form, see [3] and references therein. The block quasi-Toeplitz matrices also appear in the solution of Quasi-Birth and Death stochastic processes [8,6,7]. In this case, one of the main computational problems is solving linear systems whose input matrices are block quasi Toeplitz matrices. In [9] the authors gave a fast algorithm for solving tridiagonal quasi-Toeplitz linear systems. Inspired by [9], and other than [19,17], in this paper, we will give a fast algorithm to the block tridiagonal quasi-Toeplitz case.
2. Solving a block tridiagonal quasi-Toeplitz linear system. A block tridiagonal quasi-Toeplitz matrix is defined to be a block tridiagonal Toeplitz matrix where there are a limited number of block row changes constrained as follows: There are at most p altered block rows among the first p block rows and at most q altered block rows among the last q block rows [?], then we consider the following n-by-n nonsingular block tridiagonal quasi-Toeplitz linear system with A, X, Y are matrices of size m × m, B Hermitian matrice of size m × m, x is the nm unknown vector and f is the nm right hand side. If B = X = Y * , the matrix N becomes a symmetric block tridiagonal Toeplitz matrix [1,2,14,24]. By exploiting the structure of such matrices, many algorithms have been proposed for the solving the related linear systems [4,5,14,19].
This paper addresses the problem of solving block tridiagonal quasi-Toeplitz linear systems (i.e. B = X and B * = Y ). Inspired by Du et al. [9], we propose a more general algorithm for such systems. If we define the block diagonal matrix D as where Z = BX −1 A and V = B * Y −1 A. Therefore, the system (2.1) is equivalent to the block tridiagonal quasi-Toeplitz linear systemÑ x =f , wheref = Df.
In the following, we first recall the block LU factorization of the matrixÑ and then introduce an efficient algorithm for solving block tridiagonal quasi-Toeplitz linear systems, which combines proper matrix decomposition with the Sherman-Morrison-Woodbury formula.
2.1. Block LU decomposition. Assume that the matrixÑ admits the following decompositioñ N =LŨ , where The matrices A i and F i are lower and upper triangular matrices, respectively, and obtained by the LU factorization. Then, solving the linear system (2.1) is equivalent to solving two triangular linear systems Thus, we can give the fast algorithm to compute the block LU decomposition for solving (2.1) as follows: where I is the identity matrix of size m × m. Then, which has same structure as the matrixÑ with Equation (2.4) can be solved by a direct solver, based on the following adaptation of the algorithm presented in [18]. Then,Ñ According to Sherman-Morrison-Woodbury formula [23,12], we obtaiñ Therefore, the solution x of (2.1) is obtained as follows: where y is a solution of y = (LDU ) −1f and y 1 is the first block of the vector y, which can be efficiently computed through Algorithm 2.

Algorithm 2
Solving LDU x = y Input: L 1 , U 1 , Σ and Λ are matrices of sizes m × m Output: Finally, our algorithm for solving block tridiagonal quasi-Toeplitz linear systems can be formulated as follows:  where γ min (C) denotes the minimum eigenvalue of a matrix C with real eigenvalues and the eigenvalues of C are sorted in the ascending order:

Remark 2.4.
Note that the parameter ρ is tiny when the matrix polynomial λ 2 B * + λA + B has an eigenvalue near the unit circle.
Theorem 2.5. Assume that the computation of y = T −1f is backward stable in the sense that it is evaluated in the formỹ = (T + δ 1 ) −1f , where the perturbation δ 1 is bounded as δ 1 = O( machine ) T (see [23]). Then according to Definition 2.3, the exact solution x of the system (2.1) satisfies wherex is the computed solution, and . denotes the 2-norm.
Proof. Let Q = I + M T 1 T −1 E 1 and recall from (2.6) that then the exact solution is According to the hypothesis, we havẽ Neglecting quadratic perturbation terms, we obtain So (2.9) and (2.10) allow us to write (2.11) As a result, the forward error satisfies the estimate To derive (2.8) we observe that the perturbation matrices δ 1 and δ 2 in (2.11) possess a structure when T −1f is computed by means of the Cholesky factorization T = RR T , where R = D 1 U and D = D T 1 D 1 is also a Cholesky factorization. We again assume that R is computed with high precision and contains no rounding errors. Then, omitting quadratic perturbation terms, we can write The inequality T ≤ Ñ ⇒ RÑ −1/2 ≤ 1 and according to Theorem 2 in [17], we have that Remark 2.6. If y = O(1) x , then Remark 2.7. If the number of blocks is large, the computation of R with the double or higher precision is not expensive in comparison with the computation cost of T −1 f.

Numerical results.
4. Numerical results. In this section, we present some numerical results to illustrate the effectiveness of the proposed algorithm and we compare this algorithm with the block LU method and the Gauss algorithm [10]. The algorithms were programmed in MATLAB 9.4.0.813654 (R2018a) and the computations are done on a Intel(R) Core(TM) i7-8550U CPU @ 1.80GHz Laptop with 8.00 GB of RAM and 1.99GHz of processor. We fixed the exact solution x * = [1, 1 . . . , 1] T . The right-hand side vector was set to be f = N x * .

Conclusion.
In this paper, we proposed a generalized method for solving systems of linear equations with block tridiagonal quasi-Toeplitz matrices. Numerical examples show that the proposed algorithm is one of the good alternatives in terms of efficiency and computational time, especially when the matrix of the linear system is very large.