Band Toeplitz preconditioners for non-symmetric real Toeplitz systems by preconditioned GMRES method

Journal of Computational and Applied Mathematics xxx (xxxx) xxx

Contents lists available at ScienceDirect

Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam

Band Toeplitz preconditioners for non-symmetric real Toeplitz systems by preconditioned GMRES method✩ ∗

D. Noutsos , G. Tachyridis Department of Mathematics, University of Ioannina, GR45110, Greece

article

info

Article history: Received 30 November 2018 Received in revised form 5 April 2019 MSC: 65F08 65F10 15B05 15A18 Keywords: Non-symmetric Toeplitz matrices Band preconditioners Krylov subspace methods

a b s t r a c t In this paper we study n × n real, non-symmetric, ill conditioned Toeplitz systems of the form Tn (f )x = b. The corresponding √ generating function is a complex valued one of the form f = f1 + if2 , where i = −1, is known a priori and has roots. We note that f1 is a 2π -periodic even function while f2 is a 2π -periodic odd one. In order to solve the above system efficiently, we use the Preconditioned Generalized Minimal Residual (PGMRES) method and the Preconditioned Conjugate Gradient method applied to the Normal Equations (PCGN). We present a specific preconditioning technique that combines elimination of the roots of f and best uniform approximation or interpolation by trigonometric polynomials. The proposed preconditioner is a band Toeplitz matrix Tn (p) generated by the trigonometric polynomial p = gq. The trigonometric polynomial g is an appropriate complex polynomial having the same zeros as f while q is the trigonometric polynomial derived by the uniform approximation or interpolation of the f function g . Using Tn (p) as a preconditioner, we achieve a good clustering of the singular values of the preconditioned matrix in a small interval around 1, as well as a good clustering of the eigenvalues in a small domain in the right half-plane of the complex plane, far from zero. Moreover, we describe on how the presented preconditioning technique can be extended to the two level Toeplitz systems. Finally, we show, by various numerical experiments, that the proposed preconditioning technique can solve a Toeplitz system of the above form in a small number of iterations. © 2019 Published by Elsevier B.V.

1. Introduction In this paper we present a preconditioning technique for the solution of real, non-symmetric Toeplitz systems of the form Tn (f )x = b, where f is a 2π -periodic complex function, defined on I = [−π, π ], known a priori. From now on we will call f the generating function of the Toeplitz system since the entries tjk of the n × n matrix Tn (f ) are generated by the Fourier coefficients of f , i.e., tjk =

1 2π

∫

π

f (x)e−i(j−k)x dx, 1 ≤ j, k ≤ n, i =

√ −1.

−π

✩ This document is a collaborative effort. ∗ Corresponding author. E-mail addresses: [email protected] (D. Noutsos), [email protected] (G. Tachyridis). https://doi.org/10.1016/j.cam.2019.04.030 0377-0427/© 2019 Published by Elsevier B.V.

Please cite this article as: D. Noutsos and G. Tachyridis, Band Toeplitz preconditioners for non-symmetric real Toeplitz systems by preconditioned GMRES method, Journal of Computational and Applied Mathematics (2019), https://doi.org/10.1016/j.cam.2019.04.030.

2

D. Noutsos and G. Tachyridis / Journal of Computational and Applied Mathematics xxx (xxxx) xxx

It is easy to observe that if f is a real and even function, the generated Toeplitz matrix will be real and symmetric. If f = ig, where g is a real and odd function, then the associated Toeplitz matrix will be real and skew-symmetric. Since each real Toeplitz matrix is written uniquely as a sum of its symmetric and skew-symmetric part, the generating function of each real and non-symmetric Toeplitz matrix will be of the form f = f1 + if2 , where f1 is a real and even function and f2 is a real and odd one. Toeplitz systems arise in a variety of applications in mathematics, such as signal processing, image processing and restoration, discretization of PDEs and systems with queueing networks [1,2]. Over the years, many appropriate preconditioners have been introduced ensuring the efficiency of the Preconditioned Conjugate Gradient (PCG) method. To be more precise, R. Chan [3] in 1991 introduced a band Toeplitz preconditioner for ill-conditioned symmetric Toeplitz systems. The generating function of the preconditioner was a trigonometric polynomial g, having the roots of f . This f technique leads to the elimination of ill-conditioning, since the function g becomes positive. Later on, many authors proposed band Toeplitz preconditioners based on the elimination of ill-conditioning and on some kind of approximation of the generating function: R. Chan and P. Tang (1994) [4], S. Serra-Capizzano (1997, 1999) [5,6], D. Noutsos and P. Vassalos (2002) [7]. A preconditioner constructed as the product of a band Toeplitz matrix that eliminates the ill-conditioning and a matrix belonging to a trigonometric algebra, which leads to the superlinear convergence, was proposed by D. Noutsos and P. Vassalos (2008) [8]. We note that preconditioners from any trigonometric matrix algebra are not able to lead to superlinear convergence all by themselves [5,9,10]. In 2016 D. Noutsos, S. Serra-Capizzano and P. Vassalos [11] proposed a preconditioner belonging to τ -algebra [12,13], for Toeplitz systems generated by functions having roots of non integer order. In many applications Block Toeplitz matrices with Toeplitz Blocks (BTTB) appear [1,2,14]. BTTB matrices are called two-level Toeplitz matrices and are symbolized as Tnm (f ), where m is the dimension of each block and n is the dimension along the blocks, i.e. Tnm (f ) ∈ Rnm×nm . The generating function of such matrices is a two variate function f = f (x, y) : [−π , π]2 → C, 2π -periodic along both variables. The entries trs of the BTTB matrix are the coefficients of the Fourier expansion of f . Symmetric and positive definite BTTB systems were studied by many researchers. Efficient band BTTB preconditioners have been proposed by S. Serra [15], M. Ng [16], as well as by D. Noutsos, S. Serra-Capizzano and P. Vassalos [17–19]. Matrix algebra preconditioners have been proposed in [20] for the well-conditioned case. For the ill-conditioned case this class of preconditioners cannot be efficient, because there have been proven negative results [9,10,21–23]. Although, preconditioners based on a combination of band BTTB matrices and matrices belonging to a matrix algebra ensure the fast convergence of the PCG method [20]. We conclude that the symmetric case of Toeplitz systems is studied completely. However, the non-symmetric case seems to have plenty of open research topics. One can find theoretical results for the spectra clustering of a non-symmetric Toeplitz matrix in [24–28]. We present the cluster definitions given by Tyrtyshnikov in [27]: Definition 1 (Tyrtyshnikov). A set Φ ⊂ R is called a general eigenvalue cluster for a Hermitian (singular value for a) matrix An , An ∈ Cn×n , if for all ϵ > 0 and Φϵ being the ϵ -extension of Φ : lim

n→∞

γ n (ϵ ) n

= 0,

where γn (ϵ ) is the number of eigenvalues (singular values) of An , that lie outside of Φϵ . Φ is called a proper cluster if γn (ϵ ) ≤ c(ϵ ), where c(ϵ ) is independent of n. The ϵ -extension of Φ is the union of Φ with all the real ϵ -balls encircling Φ ’s points. For the completeness of the paper we present Theorem 4.5 of [25]: Theorem 1 (Tilli). Suppose f ∈ L2 (Q , Ch×k ) and let {Tn } be the set of block Toeplitz matrices generated by f . Then the interval [σmin (f ), σmax (f )] is a cluster for the singular values of {Tn }. We note that Q = (−π, π ) and f ∈ L2 (Q , Ch×k ) means that f is a matrix function in Ch×k with any entry fij ∈ L2 (−π, π ). In this paper we introduce a preconditioning technique for real non-symmetric Toeplitz systems, combining the idea of elimination of ill-conditioning, dividing f by a trigonometric polynomial g having the same roots with the same f order, and of the approximation of the function g by trigonometric polynomials. The elimination of ill-conditioning was inspired by [3], but here we had to deal with a complex function. We apply best uniform approximation on the f real and imaginary parts of g by trigonometric polynomials using the Remez exchange algorithm, or approximation by

trigonometric interpolation. The set of sample points consists of the Chebyshev nodes mapped on the interval [0, π ]. The paper is organized as follows: In Section 2 we analyze the spectral properties of the preconditioned Toeplitz system, using band Toeplitz preconditioners. We present a theorem for the clustering of the singular values, as well as a theorem for the location of the eigenvalues. In Section 3, we analyze the construction of the proposed preconditioner and we give bounds for the clustering of the singular values. Finally, in Section 4 we present various numerical results that show the efficiency of the proposed preconditioning technique. The last example that we present concerns the two-level case of Toeplitz systems. Please cite this article as: D. Noutsos and G. Tachyridis, Band Toeplitz preconditioners for non-symmetric real Toeplitz systems by preconditioned GMRES method, Journal of Computational and Applied Mathematics (2019), https://doi.org/10.1016/j.cam.2019.04.030.

D. Noutsos and G. Tachyridis / Journal of Computational and Applied Mathematics xxx (xxxx) xxx

3

2. Spectral properties of the preconditioned Toeplitz system — convergence of the PGMRES method We will state and prove the theoretical results concerning the rate of convergence of the PGMRES method [29,30] for the solution of real and non-symmetric Toeplitz systems, using band Toeplitz preconditioners. Our main goal is to choose an appropriate band Toeplitz preconditioner Tn (p), generated by a trigonometric polynomial p for the efficient solution of the Toeplitz system Tn (f )x = b, by the PGMRES method. To achieve this we have to choose a preconditioner such that the eigenvalues of the preconditioned matrix Tn−1 (p)Tn (f ), should belong to a domain of the open right half plane of the complex plane, far from the imaginary axis. It is well known that the rate of convergence of the PGMRES method depends on the location of eigenvalues in a small domain of the right half plane, while the convergence rate of PCGN depends on the singular value clustering in a small interval of the form [α, β ] with α < β [31]. This means ⏐ ⏐ that the trigonometric polynomial p = p1 + ip2 should be chosen so that the real part Re

( ) f p

⏐ ⏐ > 0 and the range of ⏐ pf ⏐, that is ⏐ ⏐ (⏐ ⏐) [ ⏐f ⏐ ⏐ f (x) ⏐ ⏐ ⏐ ⏐, max ⏐ range ⏐ ⏐ = min ⏐ ⏐ p

−π ≤x≤π

⏐] ⏐ ⏐ f (x) ⏐ ⏐ , ⏐ p(x) −π ≤x≤π ⏐ p(x) ⏐

should be a positive interval far from zero. The choice of the polynomial p will be extensively described in the next section concerning the construction of the preconditioner. Here we suppose that we have found such a polynomial p and we will study the spectral properties concerning the singular values as well as the eigenvalues of the preconditioned matrix. For this we state and prove the following theorems: Theorem 2. Let f ∈ L2 (−π, π ) and p be a trigonometric polynomial such that Re

⏐ ⏐ ⏐ f (x) ⏐

( ) f p

⏐ ⏐ ⏐ f (x) ⏐ > 0 and 0 < α = ess inf−π≤x≤π ⏐ p(x) ⏐≤

ess sup−π≤x≤π ⏐ p(x) ⏐ = β < ∞. Then the interval [α, β ] is a general cluster for the singular values of the preconditioned matrix Tn−1 (p)Tn (f ). Proof. It is well known that the singular values of a non-symmetric matrix A are the square roots of the eigenvalues of AT A. Thus, we will study the behavior of the eigenvalues of the matrix associated with the normal equations of the preconditioned matrix: Tn−1 (p)Tn (f )

(

)T

Tn−1 (p)Tn (f ) = TnT (f )Tn−T (p)Tn−1 (p)Tn (f )

= Tn (f¯ )Tn−1 (p¯ )Tn−1 (p)Tn (f ). The last matrix is similar to An (f , p) = Tn−1 (p)Tn (f )Tn (f¯ )Tn−1 (p¯ ). For simplicity we will study the behavior of the eigenvalues of An (f , p). Consequently, we will use the relation

( ) Tn (f ) = Tn (p)Tn

f

p

+ En ,

(1)

where En is a matrix of rank ( d)− 1, d is the bandwidth of Tn (p). Since Tn (p) is a band matrix, it is easy to observe that the f 1 matrices Tn (f ) and Tn (p)Tn p differ only in the d− first and last rows and (1) is proven. Similarly, it holds that 2 Tn (f¯ ) = Tn

(¯) f

p¯

Tn (p¯ ) + EnT .

We now study the spectrum of An (f , p). An (f , p) = Tn−1 (p)Tn (f )Tn (f¯ )Tn−1 (p¯ )

( ( ) )( (¯) ) f f = Tn−1 (p) Tn (p)Tn + En Tn Tn (p¯ ) + EnT Tn−1 (p¯ ) p p¯ ( ) (¯) (¯) ( ) f f f f −1 = Tn Tn + Tn (p)En Tn + Tn EnT Tn−1 (p¯ ) p p¯ p¯ p

(2)

+ Tn−1 (p)En EnT Tn−1 (¯p) ( ) (¯) f f = Tn Tn + Rn p p¯ where Rn has rank at most 2d − 2. Please cite this article as: D. Noutsos and G. Tachyridis, Band Toeplitz preconditioners for non-symmetric real Toeplitz systems by preconditioned GMRES method, Journal of Computational and Applied Mathematics (2019), https://doi.org/10.1016/j.cam.2019.04.030.

4

D. Noutsos and G. Tachyridis / Journal of Computational and Applied Mathematics xxx (xxxx) xxx

We recall now Theorem 4.5 of P. Tilli [25], to obtain that the interval [α, β] is a cluster of the singular values of Tn or equivalently of the square roots of the eigenvalues of Tn 4.5 of [25] is that o(n) singular values of Tn

( ) f p

( ) f p

Tn

( ) f¯ p¯

( ) f p

. The meaning of ‘‘[α, β] is a cluster’’ of Theorem

are smaller than α and all of them are smaller than β . From (2) we have

that

An (f , p) = Tn

(¯)

( ) f

Tn

p

f

p¯

+ Rn .

Thus, o(n) singular values of Tn−1 (p)Tn (f ) are less than α and at most 2d − 2 additional singular values lie outside of [α, β], which means that a constant number of singular values may be greater than β . This type of clustering has been defined by Tyrtyshnikov [27] as general cluster. Thus, we obtained the general cluster of the singular values of the preconditioned matrix Tn−1 (p)Tn (f ) and the proof is complete. □ We have to remark that Theorem 4.5 of [25] is given in a more general setting, for block Toeplitz matrices with f being matrix functions. Our settings in this paper are simpler and all the assumptions of Theorem 4.5 hold true. Let f ∈ L1 ([−π, π]) and p be a trigonometric polynomial such that ess inf−π ≤x≤π Re

Theorem 3.

(

f (x) p(x)

)

> 0. Then

the eigenvalues of ( the) preconditioned matrix Tn (p)T ( n (f)) are located in the rectangle ( R )= [α, β ] × [−γ , γ ], where α = −1

ess inf−π ≤x≤π Re

f (x) p(x)

> 0, β = ess sup−π ≤x≤π Re

f (x) p(x)

and γ = ess sup−π≤x≤π Im

f (x) p(x)

. At most 2d − 2 eigenvalues may

be outliers, which constitutes a proper cluster, where d is the bandwidth of Tn (p). The outliers belong to an O( 1n )-extension of R. Proof. We will study the range of the preconditioned matrix Tn (p)−1 Tn (f ): A=

xH Tn (p)−1 Tn (f )x

1 xH Tn (p)−1 Tn (f ) + Tn (f¯ )Tn (p¯ )−1 x

(

=

)

H x(H x 2 )x x H −1 − 1 ¯ x 1 x Tn (p) Tn (f ) − Tn (f )Tn (p¯ ) + i 2 ( ixH x ) 1 yH Tn (f )Tn (p¯ ) + Tn (p)Tn (f¯ ) y

=

Tn (f )Tn (p¯ ) − Tn (p)Tn (f¯ ) y

( H

+

(3)

yH Tn (p)Tn (p¯ )y

2 1 y i 2

)

iyH Tn (p)Tn (p¯ )y ¯ 1 y Tn (f p + pf¯ )y + yH R1 y 1 yH Tn (f p¯ − pf¯ )y + yH R3 y = + i , 2 H H 2 y Tn (|p| )y − y R2 y 2 iyH Tn (|p|2 )y − yH R2 y H

where y = Tn (p¯ )−1 x and R1 , R2 and R3 are low rank matrices of rank d − 1 each one. These matrices have non zero entries 1 at the d− first and last rows and columns, because of the band structure of Tn (p) and Tn (p¯ ). ( 2 ) It is obvious that the matrix Tn (f )Tn (p¯ ) + Tn (p)Tn (f¯ ) is symmetric while the matrix 1i Tn (f )Tn (p¯ ) − Tn (p)Tn (f¯ ) is 1 ¯ ¯ hermitian. Deleting the d− first and last rows and columns of Tn (f )Tn (p¯ ) + Tn (p)T 2 ( n (f ) as well as of Tn (f) p¯ + pf1 ) we obtain 1 ¯ that the remaining matrices coincide. The same holds true for the matrices i Tn (f )Tn (p¯ ) − Tn (p)Tn (f ) and i Tn (f p¯ − pf¯ ) as well as for the matrices Tn (p)Tn (p¯ ) and Tn (|p|2 ). Now we take the relation (3) where x is chosen from a subspace V of Rn , with dim V = n − (d − 1) and such that 1 y = Tn (p¯ )−1 x has zeros at the first and last d− entries. Then, we have: 2

( ( ( ))) f = range Re Tn p x∈V x∈V ( ( ( ))) ( ( )) f f ⊂ range Re Tn = ER Re ,

range Re Tn (p)−1 Tn (f )

(

(

))

p

x∈Rn

p

where ER(h) denotes the essential range of the function h.

(

This means that at most d − 1 eigenvalues of the real part of Tn (p)−1 Tn (f ) may lie outside the ER Re

( )) f p

, which is

the interval:

[ ( ) ( )] f (x) f (x) [α, β] = ess inf Re , ess sup Re , where α > 0. −π≤x≤π

p(x)

−π≤x≤π

p(x)

Please cite this article as: D. Noutsos and G. Tachyridis, Band Toeplitz preconditioners for non-symmetric real Toeplitz systems by preconditioned GMRES method, Journal of Computational and Applied Mathematics (2019), https://doi.org/10.1016/j.cam.2019.04.030.

D. Noutsos and G. Tachyridis / Journal of Computational and Applied Mathematics xxx (xxxx) xxx

(

Similarly, for the imaginary part we obtain that at most d − 1 eigenvalues may lie outside the ER Im

5

( )) f p

, which is the

interval:

[ ( ) ( )] f (x) f (x) [−γ , γ ] = −ess sup Im , ess sup Im . p(x)

−π ≤x≤π

p(x)

−π ≤x≤π

Combining the ranges of the real and imaginary parts we conclude that at most 2(d − 1) eigenvalues may lie outside the rectangle R = [α, β ] × [−γ , γ ], proving that it constitutes a proper cluster. The question that arises is: How far from the boundary of R, do the eigenvalues lie? In order to answer this question xH Tn (p)−1 Tn (f )x we have to estimate the Rayleigh quotient , considering that x is an eigenvector of Tn (p)−1 Tn (f ) corresponding xH x to λ. Thus,

λ=

xH Tn (p)−1 Tn (f )x

y=Tn (p¯ )−1 x

=

xH x

yH Tn (f )Tn (p¯ )y yH Tn (p)Tn (p¯ )y

.

It is known that even if Tn (f )Tn (p¯ ) and Tn (p)Tn (p¯ ) are not Toeplitz matrices, as n → ∞ they behave as the Toeplitz operators generating by f p¯ and pp¯ , respectively. To each( x ∈ [−π, )π] there corresponds an eigenvalue f (x)p¯ (x) or T p(x)p¯ (x) with the associated infinite eigenvector y = √1n 1 eix ei2x · · · [6,32,33]. Thus, for large enough n, y is close to

√1

n

(

1 eix ei2x · · · ei(n−1)x

)T

.

( ) ]T d−1 |yT2 |yT3 , where y1 , y3 ∈ R 2 and y2 ∈ Rn−(d−1) . It is obvious that ∥y1 ∥ = O √1n and ( ) ( ) ( ) ( ) ∥y3 ∥ = O √1n . We recall now relation (3) for this y. Then, obviously yH R1 y = O 1n , yH R2 y = O 1n and yH R3 y = O 1n . ( ) 1 −1 We partition y as y =

[

yT1

This means that each eigenvalue of Tn (p)

Tn (f ) belongs to an O

n

-extension of the rectangle mentioned above. □

Theorem 2 gives us that the PCGN method, because of the singular values clustering, works efficiently for the solution of the system Tn (f )x = b, while Theorem 3 guarantees the efficiency of the PGMRES method, because of the eigenvalues clustering. This is confirmed by a variety of numerical experiments shown in Section 4. Strong numerical evidence shows us that only few singular values seem to tend to zero and that the main body of outliers seems to be bounded. Also, the eigenvalues seem to be located in the rectangle described in Theorem 3, except for view ones. 3. Construction of the preconditioner As we already mentioned f = f1 + if2 is a complex function, with f1 being a real, 2π -periodic, even function and f2 a real, 2π -periodic odd one. Both f1 and f2 are defined on [−π, π]. In this section we will describe the way we construct band Toeplitz preconditioners in order to solve efficiently the Toeplitz system Tn (f )x = b by Krylov subspace methods (PGMRES, PCGN). It is obvious that the generating function of the band Toeplitz matrix is a trigonometric polynomial of the form p = p1 + ip2 , where p1 is even of degree d1 and p2 is odd of degree d2 . We will study two different cases, the well-conditioned and the ill-conditioned ones. 3.1. Well-conditioned case Well-conditioning is characterized from the fact that the function f1 is positive and |f | is bounded. This means that the range of the eigenvalues of Tn (f ) will be in a compact set in the right half plane. In this case, the band Toeplitz preconditioner will be constructed by considering some kind of approximation of the functions f1 and f2 . We propose two kinds of approximation: 1. The best uniform approximation of f1 by an even trigonometric polynomial and of f2 by an odd one, on the interval [−π , π ]. 2. Analogous trigonometric interpolation of f1 and f2 on [−π, π ]. For the best uniform approximation we apply the Remez exchange algorithm taking the nodes in a grid of k Chebyshev points of the first kind, mapped on the interval [0, π ], i.e. xj =

π 2

(

( cos

2(k − j) − 1 2k

π

)

) + 1 , j = 1, 2, . . . , k,

where k ≫ d1 , d2 . The reason we choose the nodes from the interval [0, π] is that f1 and f2 are even and odd functions, respectively on [−π , π]. We mention here that the Remez exchange algorithm does not add further cost in the PGMRES or PCGN algorithms, in order of magnitude, since it depends on k which is a constant number independent of the dimension n. In the examples that follow we choose k = 64. The following theorem guarantees the singular value clustering of the preconditioned matrix, coming out of the uniform approximation of f1 and f2 using the Remez algorithm. Please cite this article as: D. Noutsos and G. Tachyridis, Band Toeplitz preconditioners for non-symmetric real Toeplitz systems by preconditioned GMRES method, Journal of Computational and Applied Mathematics (2019), https://doi.org/10.1016/j.cam.2019.04.030.

6

D. Noutsos and G. Tachyridis / Journal of Computational and Applied Mathematics xxx (xxxx) xxx

Theorem 4. Let f = f1 + if2 , be the generating function of Tn (f ), where f1 and f2 are continuous, 2π -periodic and bounded functions. Moreover, let f1 be a positive and even function, while f2 be an odd one. Let also p = p1 + ip2 , where p1 is the best uniform even trigonometric polynomial approximation of f1 , with maximum error ϵ1 and p2 is the best uniform odd trigonometric polynomial approximation of f2 , with maximum error ϵ2 . Then, for a prescribed ϵ > 0, there exist p1 and p2 of appropriate degrees, such that the singular values of the preconditioned matrix Tn−1 (p)Tn (f ) have a general cluster in the interval [1 − M ϵ, 1 + M ϵ], where M = max−π≤x≤π

(

1 p21 (x)+p22 (x)

) 12

.

Proof. From Theorem 2 we have that the singular values of the preconditioned matrix Tn−1 (p)Tn (f ) have a general cluster ⏐ ⏐

⏐f ⏐

in the range of ⏐ p ⏐. The construction of p by the best uniform approximation gives us that f1 (x) = p1 (x) + e1 (x) and f2 (x) = p2 (x) + e2 (x),

where ∥e1 ∥∞ = ϵ1 and ∥e2 ∥∞ = ϵ2 , denote the associated errors of the best uniform approximation. Thus, ∀x ∈ [−π, π] there holds f p

=

f1 + if2 p1 + ip2

=

p1 + e1 + i(p2 + e2 ) p1 + ip2

.

T T For each x ∈ ⏐ ⏐2 e = (e1 e2 ) , as well as the vector of maxima ϵˆ = (ϵ1 ϵ2 ) . To study the ⏐ [−π ⏐ , π], we define the vector ⏐f ⏐ ⏐f ⏐ behavior of ⏐ p ⏐ we consider the function ⏐ p ⏐ :

⏐ ⏐2 2 2 2 2 ⏐f ⏐ ⏐ ⏐ = (p1 + e1 ) + (p2 + e2 ) = 1 + 2 p1 e1 + p2 e2 + e1 + e2 . ⏐p⏐ p21 + p22 p21 + p22 p21 + p22 To get lower and upper bounds we use triangle and Cauchy–Schwarz inequalities as follows:

⏐ ⏐ ⏐2 ⏐ 2 2 ⏐ ⏐f ⏐ ⏐ ⏐ ⏐ ≥ 1 − 2 ⏐ p1 e1 + p2 e2 ⏐ + e1 + e2 ⏐p⏐ ⏐ p2 + p2 ⏐ p2 + p2 2 2 1 1 e21 + e22 |p1 ||e1 | + |p2 ||e2 | + p21 + p22 p21 + p22 √ p21 + p22 √ e2 + e22 ≥1−2 2 e21 + e22 + 21 2 p1 + p2 p1 + p22 ⎛ ⎞2

≥1−2

= ⎝1 − √

1 p21 + p22

(4)

∥e∥2 ⎠ .

Thus,

⏐ ⏐ ⏐f ⏐ 1 ⏐ ⏐≥1− √ 1 ∥e∥2 ≥ 1 − max √ ∥ˆϵ ∥2 = 1 − M ∥ˆϵ ∥2 . ⏐p⏐ −π ≤ x ≤π p21 + p22 p21 + p22 For the upper bound, the reverse inequalities of (4) hold true, using + instead of − in the second term to get:

⎞2 ⎛ ⏐ ⏐2 ⏐f ⏐ 1 ⏐ ⏐ ≤ ⎝1 + √ ∥e∥2 ⎠ , ⏐p⏐ p21 + p22 and finally,

⏐ ⏐ ⏐f ⏐ 1 ⏐ ⏐≤1+ √ 1 ∥e∥2 ≤ 1 + max √ ∥ˆϵ ∥2 = 1 + M ∥ˆϵ ∥2 . ⏐p⏐ −π ≤ x ≤π p21 + p22 p21 + p22 We note that p1 should be positive for an appropriate choice of the degree d1 , since f1 is positive. This means that

( M = max

−π ≤x≤π

1 p21 (x) + p22 (x)

) 12

< ∞.

Please cite this article as: D. Noutsos and G. Tachyridis, Band Toeplitz preconditioners for non-symmetric real Toeplitz systems by preconditioned GMRES method, Journal of Computational and Applied Mathematics (2019), https://doi.org/10.1016/j.cam.2019.04.030.

D. Noutsos and G. Tachyridis / Journal of Computational and Applied Mathematics xxx (xxxx) xxx

7

The errors ϵ1 and ϵ2 of the best uniform approximation decrease as the degrees of the polynomials increase. We can choose d1 and d2 such that ∥ˆϵ ∥2 ≤ ϵ . This means that,

⏐ ⏐ ⏐f ⏐ 1 − M ϵ ≤ ⏐⏐ ⏐⏐ ≤ 1 + M ϵ, p

(5)

which is a cluster around 1. □

⏐ ⏐ ⏐f ⏐

We have to remark here that inequality (5) gives that the shape of ⏐ p ⏐ belongs to the disk with center the point (1, 0) and radius M ϵ . Theorem 3 describes well the location of the eigenvalues (rectangle). By analogous manipulations of the proof above, we obtain that this rectangle is contained in the disk with center (1, 0) and radius M ϵ . In Theorem 4 we supposed that f is continuous. In the more general setting where f is continuous in (−π, π ) but has a discontinuity at π , which means that f2 (π ) ̸ = 0, we slightly change the procedure of the approximation of f2 : We do not use the interval [0, π], but a subinterval [0, c ]. This is because any odd trigonometric polynomial has a root at the point π , while f2 may have not, which leads to a big maximum error of the approximation. Using the interval [0, c ], c < π we get a small approximation error in [0, c ] and a bigger one in the region of π . Numerical evidence shows that this does not affect the clustering of the singular values and eigenvalues. The choice of c is taken heuristically. In our numerical experiments we have chosen c = 57π . For the interpolation of f1 by an even trigonometric polynomial of degree d1 , we use the corresponding d1 +1 Chebyshev points on [0, π]. For the interpolation of f2 by an odd trigonometric polynomial of degree d2 , we use the point 0 and d2 Chebyshev points on [0, c]. Theorem 5. Let f = f1 + if2 , where f1 > 0 is a continuous, 2π -periodic function and f2 ∈ C ((−π, π )), 2π -periodic function with f2 (π ) ̸ = 0. Let also p = p1 + ip2 , where p1 is the best uniform trigonometric polynomial approximation of f1 on [−π , π], with maximum approximation error ϵ1 , and p2 is the best uniform trigonometric polynomial approximation of f2 on [−c , c ] ⊂ (−π , π ), with maximum approximation error ϵ2 , while ϵ2′ = max−π ≤x≤π |f2 (x) − p2 (x)|. Then, for a prescribed ϵ > 0, there exist p1 , p2 of appropriate degrees, such that for large enough n, L singular values, with L ≥ πc n, belong to the interval Iϵ = [1 − M ϵ, 1 + M ϵ], where M = max−π ≤x≤π

(

1 p21 (x)+p22 (x)

have a general cluster in Iϵ′ = [1 − M ϵ, 1 + M ϵ ′ ], where ϵ ′ =

√

) 12

. The remaining singular values belong outside Iϵ and

ϵ12 + ϵ2′ 2 .

Proof. It is easy to observe that all the steps of the proof of Theorem 4 hold true on the subinterval ⏐ ⏐ [−c , c ], since the

⏐f ⏐

Remez algorithm works well on this interval, to get finally the analogous clustering of the shape of ⏐ p ⏐:

⏐ ⏐ ⏐f ⏐ ⏐ ⏐ ∈ [1 − M ϵ, 1 + M ϵ], x ∈ [−c , c ]. ⏐p⏐ ⏐ ⏐ ⏐ f (x) ⏐ Let D = {x ∈ [−π , π] : ⏐ p(x) ⏐ ∈ [1 − M ϵ, 1 + M ϵ]}. Obviously, [−c , c ] ⊂ D ⊂ [−π, π]. We recall now the well-known Szegő’s theorem of the equally distribution of the singular values [24,25,27,34,35]. We consider the continuous function 0 ≤ Fh ≤ 1 as follows: Fh (z) = 1, z ≤ 1 − M ϵ − h, z ≥ 1 + M ϵ + h. Fh (z) = 0, z ∈ [1 − M ϵ, 1 + M ϵ], where h is a small positive number. Then, lim sup n→∞

1 n

#{σj ≤ 1 − M ϵ − h ∨ σj ≥ 1 + M ϵ + h}

n 1∑

π

⏐) (⏐ ⏐ f (x) ⏐ ⏐ ⏐ dx ≤ lim Fh (σj ) = Fh ⏐ n→∞ n 2π −π p(x) ⏐ j=1 ⏐) ⏐) ∫ −c (⏐ ∫ π (⏐ ⏐ f (x) ⏐ ⏐ f (x) ⏐ 1 ⏐ dx + 1 ⏐ ⏐ dx = Fh ⏐⏐ F h ⏐ 2π −π p(x) ⏐ 2π c p(x) ⏐ (∫ −c ) ∫ π 1 2(π − c) π −c 1dx + 1dx = = , ≤ 2π 2π π −π c 1

∫

where # denotes the cardinality of a set and ∨ the logical OR. The last inequality is due to the fact that [−c , c ] ⊂ D. Taking h > 0 tending to zero we get that lim sup n→∞

1 n

#{σj ≤ 1 − M ϵ ∨ σj ≥ 1 + M ϵ} ≤

π −c . π

Please cite this article as: D. Noutsos and G. Tachyridis, Band Toeplitz preconditioners for non-symmetric real Toeplitz systems by preconditioned GMRES method, Journal of Computational and Applied Mathematics (2019), https://doi.org/10.1016/j.cam.2019.04.030.

8

D. Noutsos and G. Tachyridis / Journal of Computational and Applied Mathematics xxx (xxxx) xxx

Thus, lim sup n→∞

1 n

#{σj ∈ [1 − M ϵ, 1 + M ϵ]} ≥

c

π

.

This means that for large enough n we get that L singular values belong to the interval Iϵ = [1 − M ϵ, 1 + M ϵ], with L ≥ πc n. We will show that the rest of the singular values have a general cluster on Iϵ ′ = [1 − M ϵ, 1 + M ϵ ′ ] and the main

√

body of them belongs to [1 + M ϵ, 1 + M ϵ ′ ], where ϵ ′ =

ϵ12 + ϵ2′ 2 .

We observe that p2 (π ) = 0 while f2 (π ) ̸ = 0. Due to the continuity arguments, we can easily conclude that |f2 | > |p2 | and sign(f2 ) = sign(p2 ) in a region of π . We consider that the value of c is chosen close enough to π , such that the relations above hold true. For the analysis we define the function s as follows:

{ s(z) =

1,z ≥ 0 −1 , z < 0

For each x ∈ [−π , π] it holds,

−ϵ1 ≤ e1 (x) ≤ ϵ1 , x ∈ [−π, π] −ϵ2 ≤ e2 (x) ≤ ϵ2 ⇔ −ϵ2 ≤ s(p2 (x))e2 (x) ≤ ϵ2 , x ∈ [−c , c ].

(6)

From the assumption above, we obtain that in the interval [c , π], the function e2 has the same sign as p2 . The same holds true in the interval [−π, −c ]. Thus, 0 ≤ s(p2 (x))e2 (x) ≤ ϵ2′ , |x| ∈ [c , π].

(7)

Combining the relations (6) and (7) we obtain that

− ϵ2 ≤ s(p2 (x))e2 (x) ≤ ϵ2′ , x ∈ [−π, π]. ⏐ ⏐2 ⏐f ⏐ We estimate now the bounds of ⏐ p ⏐ :

(8)

⏐ ⏐2 2 2 2 2 ⏐f ⏐ ⏐ ⏐ = (p1 + e1 ) + (p2 + e2 ) = (p1 + e1 ) + (s(p2 )p2 + s(p2 )e2 ) ⏐p⏐ p2 + p2 p2 + p2 2

1

=

2

1

(p1 + e1 )2 + (|p2 | + s(p2 )e2 )2 p21 + p22

.

As in the proof of Theorem 4 we define the vectors e = (e1 e2 )T , ϵˆ = (ϵ1 ϵ2 )T and ϵˆ ′ = (ϵ1 ϵ2′ )T . From (6), (8) and the fact that p1 > 0 we get: (p1 − ϵ1 )2 + (|p2 | − ϵ2 )2 p21 + p22

⏐2 ′ 2 2 ⏐ ⏐ ≤ (p1 + ϵ1 ) + (|p2 | + ϵ2 ) . ⏐ 2 2 p p +p

⏐ ⏐f ≤ ⏐⏐

2

1

The lower bound gives: (p1 − ϵ1 )2 + (|p2 | − ϵ2 )2 p21 + p22

=1−2

p1 ϵ1 + |p2 |ϵ2 p21 + p22

+

ϵ12 + ϵ22 p21 + p22

√ ≥1−2

p21 + p22

p21

+

p22

∥ˆϵ ∥2 +

∥ˆϵ ∥22 ≥ (1 − M ∥ˆϵ ∥2 )2 , + p22

p21

while the upper bound gives: (p1 + ϵ1 )2 + (|p2 | + ϵ2′ )2 p21 + p22

≤ (1 + M ∥ˆϵ ′ ∥2 )2 .

Thus,

⏐ ⏐ ⏐f ⏐ 1 − M ϵ ≤ 1 − M ∥ˆϵ ∥2 ≤ ⏐⏐ ⏐⏐ ≤ 1 + M ∥ˆϵ ′ ∥2 = 1 + M ϵ ′ , p and the proof is completed.

□

We have to remark here that the essential difference of the two intervals above is that ϵ tends to zero as the degrees of the polynomials increase, while ϵ ′ tends to a constant greater than zero. In practice, ϵ becomes bigger as c is chosen closer to π . This means that we get fewer outliers from the interval [1 − M ϵ, 1 + M ϵ], but this interval becomes bigger. That is the reason of choosing c heuristically. Please cite this article as: D. Noutsos and G. Tachyridis, Band Toeplitz preconditioners for non-symmetric real Toeplitz systems by preconditioned GMRES method, Journal of Computational and Applied Mathematics (2019), https://doi.org/10.1016/j.cam.2019.04.030.

D. Noutsos and G. Tachyridis / Journal of Computational and Applied Mathematics xxx (xxxx) xxx

9

3.2. Ill-conditioned case In contrast to the well conditioned case, in the ill-conditioned one the function f1 has roots. Since f2 , as an odd function, has in any case a root at zero, we first study the case where f1 has also a root at zero. 3.2.1. f has a single root at zero We remark that f1 is an even function and in the current case f1 does not change sign on [−π, π]. Suppose that m1 and m2 are the multiplicities of the root of f1 and f2 , respectively. We will study the case where m1 is an even integer, while m2 is an odd one. Suppose first that m1 < m2 , which means that the multiplicity of the root of the function f is m1 . Then, following the well-known technique proposed by Raymond Chan in [3] we eliminate the ill-conditioning, dividing by the trigonometric polynomial g = (2 − 2 cos (x))(

m1 ) 2

.

Thus, we obtain the function fˆ =

f g

=

f1 + if2

f1

=

g

g

+i

f2 g

= fˆ1 + ifˆ2 .

Now, we get fˆ1 > 0 (the well conditioned case). Hence, we can use the same approximation techniques for fˆ instead of f . Let q = q1 + iq2 be the associated approximation of fˆ . Dividing by this polynomial we get fˆ q

=

f g

q

=

f1 +if2 g

q1 + iq2

=

f1 + if2 gq1 + igq2

=

f p

. f

⏐ ⏐p constitutes ⏐f ⏐ −1 a cluster around one and from Theorem 2 the singular values of Tn (p)Tn (f ) are clustered in the range of ⏐ p ⏐ in the sense We now choose the band Toeplitz matrix Tn (p) as a preconditioner, where p = gq1 + igq2 . The range of

of a general cluster. We have to remark here that if f is a continuous function, Theorem 4 holds true for the preconditioned matrix Tn−1 (p)Tn (f ). On the other hand, if f2 has a discontinuity at π , Theorem 5 holds true. We remark also that fˆ2 has still a root at zero, but this does not have any significant effect since fˆ2 is an odd function. In the case where m2 < m1 , the multiplicity of f is m2 , which corresponds to the imaginary direction. If we try to m1 eliminate the ill-conditioning, using the same technique as above, that is dividing by (2 − 2 cos (x))( 2 ) , the imaginary term will tend to infinity at zero. On the other hand, if we try to eliminate the effect of the root of f2 dividing by i (sin (x))m2 , then the imaginary term of f becomes real one without roots and the real one becomes imaginary. The problem in this case is that the becoming real part tends to infinity as x → ±π . To avoid this we divide by a combination of both functions and more precisely by g = g1 + ig2 = (2 − 2 cos (x))

m1 2

+ i (sin (x))m2 .

Then, we get the function fˆ =

f1 + if2 g1 + ig2

=

f1 g1 + f2 g2 g12 + g22

+i

f2 g1 − f1 g2 g12 + g22

= fˆ1 + ifˆ2 .

Now the real part fˆ1 is positive and bounded while the imaginary one fˆ2 has a root at zero. Thus we lead to the well conditioned case. Following the same procedure for the approximation of fˆ as in the previous case, we construct the polynomial q = q1 + iq2 . Dividing by this polynomial we get fˆ q

=

f g

q

=

f gq

=

f p

,

where p = gq = (g1 + ig2 )(q1 + iq2 ) = g1 q1 − g2 q2 + i(g1 q2 + g2 q1 ) = p1 + ip2 . As we described previously, the clustering of the singular values is ensured and the efficiency of the PCGN method is guaranteed. It is easy to check that Theorems 4 and 5 hold true, under the associated assumptions. Theorem 3 guarantees the efficiency of the PGMRES method. The assumption that m1 is an even integer and m2 is an odd one is necessary. There are even functions having a root of odd multiplicity at zero (e.g. f1 = |x|), but their derivative is not defined at zero. The same holds true for f2 (e.g. f2 = |x|x). The smoothness is required at the point of the root. This case cannot be covered and an efficient band Toeplitz preconditioner does not exist. Please cite this article as: D. Noutsos and G. Tachyridis, Band Toeplitz preconditioners for non-symmetric real Toeplitz systems by preconditioned GMRES method, Journal of Computational and Applied Mathematics (2019), https://doi.org/10.1016/j.cam.2019.04.030.

10

D. Noutsos and G. Tachyridis / Journal of Computational and Applied Mathematics xxx (xxxx) xxx

3.2.2. f has a root at a point different from zero Suppose that f has a root at the point x0 ∈ [−π, π ), x0 ̸ = 0. Thus, x0 is a root of f1 and f2 with multiplicity m1 and m2 , respectively. Here, it is not necessary that m1 is an even number and m2 is an odd one. Since f1 is an even function and f2 an odd one, −x0 is also a root of both functions with the multiplicities as in x0 . The function f2 has an additional root at zero with odd multiplicity. Since ±x0 are roots of f1 of order m1 , in small regions of ±x0 : Iϵ = [−ϵ − x0 , −x0 + ϵ] ∪ [x0 − ϵ, x0 + ϵ], the function f1 should have the form f1 (x) = c1 (x)(x − x0 )m1 (x + x0 )m1 + o (x − x0 )m1 (x + x0 )m1 ,

)

(

(9)

where c1 (x) is a bounded function far from zero which conserves the sign in Iϵ . We have to remark here that the function f1 should be smooth enough at the points of zeros. Especially it should be smooth of order m1 . For example, the even function |x − 1||x + 1| seems to have roots of multiplicity equal to one at the points ±1, but it is not differentiable at those points. This function cannot be written in the form (9) and its roots cannot be eliminated. However, when f1 is smooth of order m1 , its roots can be eliminated dividing by the trigonometric polynomial

(

(

sign(c1 (x)) sin 1

= sign(c1 (x)) 2 The coefficient

1

m1

m1 2

x − x0 2

))m1 ( sin

(

x + x0

))m1

2

(cos (x0 ) − cos (x))m1 .

does not play any role thus we can use the simplest form:

2 2

sign(c1 (x)) (cos (x0 ) − cos (x))m1 .

(10)

We have to remark here that we could not use the trigonometric polynomial sign(c1 (x)) (sin (x − x0 ))m1 (sin (x + x0 ))m1 , since it has two additional roots |x0 | − π and −|x0 | + π in (−π, π ). The same analysis holds true for the function f2 , but we have an additional odd root at zero of multiplicity m0 . Finally, the appropriate trigonometric polynomial for the elimination of the roots of f2 is q2 (x) = sign(c2 (x)) (cos (x0 ) − cos (x))m2 (sin (x))m0 ,

(11)

where c2 (x) plays the same role as c1 (x) in (9) but in the set Iϵ′ = [−ϵ − x0 , −x0 + ϵ] ∪ [−ϵ, ϵ] ∪ [x0 − ϵ, x0 + ϵ]. The smoothness of the function f2 is also required as it was described above, for f1 . In the case where m1 ≤ m2 , we need to eliminate the roots of f1 . Thus g is the trigonometric polynomial given by (10). If m1 < m2 the roots of f2 at ±x0 are remaining as roots of fewer multiplicity. Thus, in this special case, the smoothness of f2 is not required. For example, if f = (x − 1)(x + 1) + ix|x − 1|(x − 1)|x + 1|(x + 1), although f2 does not have the required smoothness, the procedure of the construction of the preconditioner, works well. The same holds true when m2 is not an integer number. In the case where m2 < m1 we have to use the combination of the trigonometric polynomials given by (10) and (11). Thus we have, g = g1 + ig2 = sign(c1 (x)) (cos (x0 ) − cos (x))m1

(12)

+ i sign(c2 (x)) (cos (x0 ) − cos (x))m2 (sin (x))m0 . f

Then we follow the same technique to approximate fˆ = g by a trigonometric polynomial and finally, the band Toeplitz preconditioner will be generated by p = gq and the theory presented above holds true for the preconditioned matrix Tn−1 (p)Tn (f ). 3.2.3. Roots of f1 and f2 at different points Suppose that f1 has roots at ±x1 ∈ [−π, π ) of order m1 and f2 has roots at ±x2 ∈ [−π, π ) of order m2 and an additional root at zero of order m0 . It is easy to observe that f does not have any root but f takes values in both half planes of the complex plane. In order to eliminate the roots of f1 and construct an efficient preconditioner we divide by a more specific form of (12), which can be seen below: g = g1 + ig2 = sign(c1 (x)) (cos (x1 ) − cos (x))m1

+ i sign(c2 (x)) (cos (x2 ) − cos (x))m2 (sin (x))m0 .

(13)

Please cite this article as: D. Noutsos and G. Tachyridis, Band Toeplitz preconditioners for non-symmetric real Toeplitz systems by preconditioned GMRES method, Journal of Computational and Applied Mathematics (2019), https://doi.org/10.1016/j.cam.2019.04.030.

D. Noutsos and G. Tachyridis / Journal of Computational and Applied Mathematics xxx (xxxx) xxx

11

Dividing by the trigonometric polynomial given by (13), we eliminate the roots of f1 and the range of fˆ1 belongs to a bounded set of the positive halfplane. We note that if we try to eliminate the roots of f1 , dividing by a trigonometric polynomial analogous to (10) i.e. sign(c1 (x)) (cos (x1 ) − cos (x))m1 , the imaginary term would tend to infinity at ±x1 . We remark that if f1 has roots at ±x1 and f2 has only the root zero (x2 = 0), we choose as g the trigonometric polynomial g = sign(c1 (x)) (cos (x1 ) − cos (x))m1 + i sign(c2 (x)) sin (x)m0 . 3.2.4. Roots of f at multiple points Suppose that the function f1 has roots at the points ±x1 , ±x2 , . . . , ±xk of order m1 , m2 , . . . , mk , respectively and f2 has also the same roots of order ℓ1 , ℓ2 , . . . , ℓk and an additional root at zero of multiplicity ℓ0 . If mi ≤ ℓi , ∀i = 1, 2, . . . , k, we can choose the trigonometric polynomial which eliminates all the roots of f1 . This will constitute a product of trigonometric polynomials given by (10) for all the roots, i.e. g = sign(c1 (x))

k ∏

(cos (xi ) − cos (x))mi .

(14)

i=1

If f1 has also a root at zero with multiplicity m0 ≤ ℓ0 , the trigonometric polynomial which eliminates all the roots of f1 has the form: g = sign(c1 (x))(2 − 2 cos (x))

m0 2

k ∏

(cos (xi ) − cos (x))mi .

i=1

The reason for which we demand m0 ≤ ℓ0 is that otherwise, the imaginary part of fˆ would tend to infinity at zero. The case where m0 > ℓ0 can be covered by a different choice of g, as it is described below. We define the set of indexes Q = {i : mi ≤ ℓi } and R = {i : mi > ℓi }. Obviously, if R = ∅, then #Q = k and this case is covered above. We will describe the case where R ̸ = ∅. In this case we cannot use the trigonometric polynomial given by (14), because at the points of roots ±xi : i ∈ R we will cause the imaginary part to tend to infinity. Thus, we have to use a combination of trigonometric polynomials in order to eliminate both the roots of f1 and f2 . It is obvious thatf has roots at the points ±xi with multiplicity min {mi , ℓi }. Thus, it is necessary that g should have the function gˆ =

k ∏

(cos (xi ) − cos (x))min {mi ,ℓi } ,

i=1

as a term in its product. The other term should eliminate the remaining roots. The remaining roots of f1 are at different points than those of f2 , since we reduced the multiplicity of the roots (of f ) by min {mi , ℓi }. The remaining roots of f1 correspond to indexes i ∈ R, while those of f2 correspond to indexes i ∈ Q and we lead to a generalization of the case of Section 3.2.3. Thus, taking into account the relation (13) this term should be: g˜ = sign(c1 (x))

∏

(cos (xi ) − cos (x))mi −ℓi

i∈R

+ i sign(c2 (x)) (sin (x))ℓ0

∏

(cos (xi ) − cos (x))ℓi −mi .

i∈Q

The product of gˆ and g˜ gives us the function g which, after some manipulations, is written as: g = sign(c1 (x))

k ∏

(cos (xi ) − cos (x))mi

i=1

ℓ0

+ i sign(c2 (x)) (sin (x))

k ∏

(15) ℓi

(cos (xi ) − cos (x)) .

i=1

We introduced the appropriate trigonometric polynomial for the case where the function f has roots at multiple points. This choice of g (i.e. (15)) covers also the case where the function f has roots at multiple points as above and there may also be some points where f1 has zeros and f2 does not and vice versa. Let xj ̸ = 0 be a point, where f1 has a root of multiplicity mj and f2 does not. We assume that ℓj = 0. Now let xj be a point, where f2 has root of multiplicity ℓj and f1 does not. We assume that mj = 0. As it was mentioned above, this case is covered by the relation (15), where some mi ’s and ℓi ’s may be equal to zero. m0 We have to remark that if f1 has a root at zero of multiplicity m0 , we multiply the first term of (15) by (2 − 2 cos (x)) 2 , which is the one that eliminates this root. After this, all the roots will be eliminated and we follow the well-known f technique of the approximation of the function fˆ = g . Please cite this article as: D. Noutsos and G. Tachyridis, Band Toeplitz preconditioners for non-symmetric real Toeplitz systems by preconditioned GMRES method, Journal of Computational and Applied Mathematics (2019), https://doi.org/10.1016/j.cam.2019.04.030.

12

D. Noutsos and G. Tachyridis / Journal of Computational and Applied Mathematics xxx (xxxx) xxx Table 1 Number of iterations when f (x) = x2 + 1 + ih1 (x). n

256 512 1024 2048

PGMRES

PCGN

In

R4,4

R6,6

R8,6

In

R4,4

R6,6

R8,6

31 30 29 29

8 8 8 8

7 7 7 6

6 6 6 6

72 74 73 72

37 31 30 30

34 30 30 29

38 31 29 30

3.3. Two-level case The proposed technique above could be extended, under suitable modifications, to the solution of non-symmetric and non-definite two-level Toeplitz systems by PGMRES method. Although, in the analysis of this case some further difficulties take place. The first one is that no best uniform approximation exists for two variate functions. We can overcome this difficulty by using interpolation in two variate functions. The other difficulty is that f (f1 and/or f2 ) may have curves of roots, instead of simple points. This case can be covered by using analogous analysis as in [17]. However, in the simple case where the function f is of separated variables all the analyses above can be straightforward generalized. That is because in this case the BTTB matrix is produced as a sum of tensor products of unilevel Toeplitz matrices. After the separation of the variables we can use the Remez algorithm for each univariate function and proceed to the construction of the preconditioner using the associated tensor products. We show the validity of the theory above by characteristic Example 6. 4. Numerical experiments In this section, we present a variety of numerical experiments in order to show the efficiency of our proposed technique of preconditioning. The experiments were carried out using Matlab. In all of them the right hand side vector b, of the system Tn (f )x = b, was chosen such that the solution of the system to be the vector of all ones (1 1 · · · 1)T . Our initial ∥r (k) ∥ guess was the zero vector and as stopping criterion we have chosen: ∥r (0) ∥2 ≤ 10−6 , where r (k) is the residual vector in 2

the kth iteration and r (0) = b. In the tables of iterations below, we use the following notation: By In we denote that no preconditioning took place, B denotes that the preconditioner is the band Toeplitz matrix Tn (g) and Rd1 ,d2 denotes that the preconditioner is Tn (p) f derived by the best uniform approximation of g . We note that p = gq with q = q1 + iq2 and d1 denotes the degree of q1 , while d2 the degree of q2 . Following the same formalism, by Ind1 ,d2 we denote that the preconditioner was constructed interpolating fˆ1 and fˆ2 by trigonometric polynomials of degrees d1 and d2 , respectively. We also present the number of the singular values which lie outside the interval [1 − M ϵ, 1 + M ϵ] (see Theorems 4, 5), as SV − out.

⎧ π ⎪ ⎪ ⎪−π − x , −π ≤ x < − ⎪ ⎪ 2 ⎨ π . It is obvious that h is a 2π -periodic and π Example 1. Let f = x2 + 1 + ih1 (x), where h1 (x) = 1 x ,− ≤ x < ⎪ 2 2 ⎪ ⎪ ⎪ ⎪ ⎩ π −x , π ≤x≤π

2 continuous function. We also note that f1 = x2 + 1 is a positive function. The number of iterations needed in order to have convergence using PGMRES and PCGN, is shown in Table 1. We observe that PGMRES is much faster than PCGN. This is explained by the type of clustering that characterizes the convergence of the two methods. We have general cluster for the singular values, some of them may be close to zero, while proper cluster for the eigenvalues in a rectangle far from the origin. Fig. 1 shows the spectra clustering, when n = 2048. The blue lines in Fig. 1a denote the values of 1 − M ϵ and 1 + M ϵ . We note that f1 was approximated by p1 of degree 8 and h1 by p2 of degree 6. As we can see all the singular values are located between 1 − M ϵ and 1 + M ϵ and Theorem 4 holds true. In Fig. 1b we can observe that the eigenvalues of the preconditioned system lie inside the rectangle [0.822, 1.178] × [−0.125, 0.125], f that is the one having edges bounded by the real and imaginary parts of p and Theorem 3 holds true. This property is characterized by the proper clustering of the eigenvalues [36]. Example 2. Let f = x2 + ix3 . It is obvious that f1 = x2 , f2 = x3 and both f1 and f2 have a root at zero with multiplicities m1 = 2 and m2 = 3, respectively. Thus, we eliminate the root of f dividing by the trigonometric polynomial g(x) = 2 − 2 cos (x), as it is described in Section 3.2.1. Then, we approximate the function construct the preconditioner.

f g

as described in Section 3 and

Please cite this article as: D. Noutsos and G. Tachyridis, Band Toeplitz preconditioners for non-symmetric real Toeplitz systems by preconditioned GMRES method, Journal of Computational and Applied Mathematics (2019), https://doi.org/10.1016/j.cam.2019.04.030.

D. Noutsos and G. Tachyridis / Journal of Computational and Applied Mathematics xxx (xxxx) xxx

13

Fig. 1. Spectra when f (x) = x2 + 1 + ih1 (x). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 2. x2 + ix3 . (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) Table 2 Number of iterations and outliers when f (x) = x2 + ix3 . n

256 512 1024 2048

PGMRES

PCGN

SV-out

In

B

R4,4

R6,6

In

B

R4,4

R6,6

256 512 1022 2029

67 70 69 68

24 27 28 28

22 26 27 27

– – – –

80 93 104 115

37 43 47 52

35 41 44 48

Fig. 2 shows the real and imaginary parts of the function f = x2 + ix3 (black line), fˆ =

58 114 226 450

f g

(blue line), q: the best uniform

trigonometric approximation using polynomials of 6th degree for fˆ1 and fˆ2 (orange line), and

f gq

(green line).

The number of iterations and outliers (with respect to 1 ± M ϵ ) is shown in Table 2. We note that the outliers concern R6,6 . We approximate fˆ2 on the interval [− 57π , 57π ] (i.e. we choose c = 57π ). We observe that the number of outliers concerning Iϵ = [1 − M ϵ, 1 + M ϵ] = [0.931, 1.069], is smaller than n − πc n and Theorem 5 holds true. For instance, when n = 2048, we have 450 outliers, while [n − πc n] = 585. The spectra clustering for n = 2048 is shown in Fig. 3: Fig. 3a shows the singular value clustering, Fig. 3b shows the last 50 singular values and Fig. 3c shows the eigenvalue clustering inside the rectangle [0.933, 2.015] × [−3.236, 3.236]. In Fig. 3a, b the orange lines denote the interval Iϵ and the blue one the value 1 + M ϵ ′ . The red stars denote the singular values which lie outside Iϵ′ = [1 − M ϵ, 1 + M ϵ ′ ] = [0.931, 8.535] and characterize the general cluster of Theorem 2. Fig. 3c shows the validity of Theorem 3. It is remarkable that the main mass of the eigenvalues is clustering so close to the point (1, 0). In the examples that follow we concentrate only on the behavior of PGMRES since we have noticed that it is much faster than PCGN. Please cite this article as: D. Noutsos and G. Tachyridis, Band Toeplitz preconditioners for non-symmetric real Toeplitz systems by preconditioned GMRES method, Journal of Computational and Applied Mathematics (2019), https://doi.org/10.1016/j.cam.2019.04.030.

14

D. Noutsos and G. Tachyridis / Journal of Computational and Applied Mathematics xxx (xxxx) xxx

Fig. 3. Spectra when f (x) = x2 + ix3 . (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 4. x2 + ix. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Example 3. Let f = x2 + ix. In this example f1 = x2 , f2 = x and both f1 and f2 have a root at zero with multiplicities m1 = 2 and m2 = 1, respectively. Thus, we will choose g(x) = 2 − 2 cos (x) + i sin (x), as it is described in Section 3.2.1. Then, we approximate/interpolate the function

f g

(by trigonometric polynomials) and we construct the preconditioner.

Fig. 4 shows the real and imaginary parts of f , fˆ =

f , g

q concerning R4,4 and

f , gq

as in the previous example.

Please cite this article as: D. Noutsos and G. Tachyridis, Band Toeplitz preconditioners for non-symmetric real Toeplitz systems by preconditioned GMRES method, Journal of Computational and Applied Mathematics (2019), https://doi.org/10.1016/j.cam.2019.04.030.

D. Noutsos and G. Tachyridis / Journal of Computational and Applied Mathematics xxx (xxxx) xxx

15

Table 3 Number of iterations and outliers when f (x) = x2 + ix. n

In

B

R4,4

In4,4

R10,10

In10,10

SV-out

256 512 1024 2048

256 512 1024 2048

11 11 10 10

6 6 6 6

6 6 6 5

6 6 5 5

12 12 12 11

2 2 2 2

Fig. 5. Spectra when f (x) = x2 + ix. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) Table 4 Number of iterations and outliers when f (x) = x2 − 1 + ih2 (x). n

In

B

R4,4

SV-out

256 512 1024 2048

256 512 1024 2047

15 15 16 15

6 6 6 6

4 4 4 4

The number of iterations is shown in Table 3. Moreover, Table 3 shows the number of singular value outliers concerning R4,4 . We observe that the preconditioners derived by best uniform approximation and interpolation have almost the same effectiveness when the degrees of the polynomials q1 and q2 are small. However, as the degrees become bigger, the preconditioner derived by best uniform approximation works better, while the one derived by interpolation is not effective enough, due to the variation of the trigonometric polynomial (arising from interpolation). Thus, it is confirmed that the best uniform approximation gives more efficient preconditioners. The spectra clustering concerning n = 2048 and R4,4 , is shown in Fig. 5: Fig. 5a shows the clustering of the singular values, while Fig. 5b the clustering of the eigenvalues. In Fig. 5a the orange lines denote the interval Iϵ = [0.904, 1.096] and the blue one the value of 1 + M ϵ ′ = 1.809. As in the previous example, with red stars we denote the outliers that form the general cluster. Fig. 5b shows the eigenvalue clustering inside the rectangle [0.915, 1.095] × [−0.331, 0.331].

⎧ 1 ⎪ ⎪ ⎪−1 − x , −π ≤ x < − ⎪ ⎪ 2 ⎨ 1 1 Example 4. Let f = x2 − 1 + ih2 (x), where h2 (x) = . The function f1 has single roots at ±1 of x ,− ≤ x < ⎪ 2 2 ⎪ ⎪ ⎪ ⎪ ⎩ 1−x , 1 ≤x≤π 2 order m1 = 1 and f2 = h2 has also roots at ±1 of order m2 = 1 and a root at zero of order m0 = 1. Thus, we choose f g = cos (1) − cos (x), as it is described in Section 3.2.2. Then, we approximate g by trigonometric polynomials of 4th degree and construct the preconditioner. The number of iterations and outliers is shown in Table 4. The spectra clustering is shown in Fig. 6. Fig. 6a shows the clustering of the singular values when n = 2048. Fig. 6b shows ( the ) clustering of the eigenvalues. We observe that two eigenvalues have real part greater than ess sup−π ≤x≤π

Re

f (x) p(x)

, characterizing the proper clustering in the rectangle [0.888, 1.088] × [−0.253, 0.253].

Please cite this article as: D. Noutsos and G. Tachyridis, Band Toeplitz preconditioners for non-symmetric real Toeplitz systems by preconditioned GMRES method, Journal of Computational and Applied Mathematics (2019), https://doi.org/10.1016/j.cam.2019.04.030.

16

D. Noutsos and G. Tachyridis / Journal of Computational and Applied Mathematics xxx (xxxx) xxx Table 5 Number of iterations when f (x) = (x2 − 1)2 + ix(x2 − 4). n

In

B

R8,6

256 512 1024 2048

151 197 223 229

25 25 25 24

12 11 11 11

Fig. 6. Spectra when f (x) = x2 − 1 + ih2 (x). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 7. Eigenvalues clustering when f (x) = (x2 − 1)2 + ix(x2 − 4). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Example 5. Let f = (x2 − 1)2 + ix(x2 − 4). Obviously, in this example the real and imaginary parts of f have roots at different points since f1 = (x2 − 1)2 and f2 = x(x2 − 4). More specifically, f1 has roots at ±1 with multiplicity m1 = 1 and f2 has roots at 0 and ±2 with multiplicities m0 = m2 = 1. Thus, we choose g = (cos (1) − cos (x))2 + i sin (x) (cos (2) − cos (x)), as it is described in Section 3.2.3. Then, we construct the preconditioner R48,6 . The number of iterations using PGMRES is shown in Table 5. The eigenvalues proper clustering inside the rectangle [0.684, 1.953] × [−0.246, 0.246] is shown in Fig. 7. Example 6. Let f = f (x, y) = x2 + y2 + i(x + y). Obviously, f1 (x, y) = x2 + y2 and f2 (x, y) = x + y. We observe that f is a function of separated variables. The function f1 has a root at the point (0, 0) with multiplicity 2 and f2 has also a root at (0, 0) with multiplicity 1. Let h(z) = z 2 + iz , z ∈ [−π, π], then the function f is separated as: f (x, y) = h(x) + h(y). We will eliminate the root of h(z), dividing by g(z) = 2 − 2 cos (z) + i sin (z). Then, we approximate the function gh by the trigonometric polynomial q = q1 + iq2 , where q1 and q2 are polynomials of 4th degree. Then, we construct the unilevel band Toeplitz matrix Tn (p) = Tn (gq), which is the corresponding preconditioner for Tn (h). Finally, the two-level band Toeplitz preconditioner is constructed by the tensor product Tnm (pˆ ) = In ⊗ Tm (p) + Tn (p) ⊗ Im , where pˆ (x, y) = p(x) + p(y). The number of iterations for different block dimensions is shown in Table 6. The spectra clustering is shown in Fig. 8. Fig. 8a shows the clustering of the singular values when n = m = 128. Fig. 8b shows the clustering of the eigenvalues. Please cite this article as: D. Noutsos and G. Tachyridis, Band Toeplitz preconditioners for non-symmetric real Toeplitz systems by preconditioned GMRES method, Journal of Computational and Applied Mathematics (2019), https://doi.org/10.1016/j.cam.2019.04.030.

D. Noutsos and G. Tachyridis / Journal of Computational and Applied Mathematics xxx (xxxx) xxx

17

Table 6 Number of iterations when f (x, y) = x2 + y2 + i(x + y). n

m 16

32

64

128

16 32 64 128

6 6 6 6

6 6 6 6

6 6 6 6

6 6 6 6

Fig. 8. Spectra when f (x, y) = x2 + y2 + i(x + y).

5. Conclusion We have studied the problem of the solution of real and non-symmetric Toeplitz systems by Krylov subspace iterative methods and especially by PGMRES and PCGN. We proposed a procedure of construction of the band Toeplitz preconditioner in both well-conditioned and ill-conditioned cases. This procedure is based on the elimination of the roots of the generating function and on the trigonometric polynomial approximation. Spectral properties have been proven and cover theoretically the efficiency of the proposed preconditioners. We note that spectrally equivalent systems arise during the discretization process of Fractional Differential Equations using Grünwald formula [37]. A variety of numerical experiments have shown the validity of the theoretical results. It is well known that in real and scientific applications, two-level Toeplitz systems mostly appeared. We strongly believe that our technique can be generalized and extended to the two-level case, as it is shown in characteristic Example 6. Nevertheless, additional difficulties are appeared. This case is an open topic and needs further study. Our presented work provides the theoretical background and the basic ideas to construct preconditioners that will help in this further study. References [1] R.H. Chan, X.-Q. Jin, An Introduction to Iterative Toeplitz Solvers, SIAM, Philadelphia, 2007, http://dx.doi.org/10.1137/1.9780898718850. [2] M. Ng, Iterative Methods for Toeplitz Systems, Oxford University Press, New York, 2004. [3] R.H. Chan, Toeplitz preconditioners for toeplitz systems with nonnegative generating functions, IMA J. Numer. Anal. 11 (1991) 333345, http://dx.doi.org/10.1093/imanum/11.3.333. [4] R.H. Chan, T. Tang, Fast band–Toeplitz preconditioners for Hermitian Toeplitz systems, SIAM J. Sci. Comput. 15 (1994) 164–171, http: //dx.doi.org/10.1137/0915011. [5] S. Serra, Optimal, quasi-optimal and superlinear band-Toeplitz preconditioners for asymptotically ill-conditioned positive definite Toeplitz systems, Math. Comp. 66 (1997) 651–665, http://dx.doi.org/10.1090/S0025-5718-97-00833-8. [6] S. Serra, How to choose the best iterative strategy for symmetric Toeplitz systems, SIAM J. Numer. Anal. 36 (1999) 1078–1103, http: //dx.doi.org/10.1137/S0036142996311866. [7] D. Noutsos, P. Vassalos, New band Toeplitz preconditioners for ill-conditioned symmetric positive definite Toeplitz systems, SIAM J. Matrix Anal. Appl. 23 (2002) 728–743, http://dx.doi.org/10.1137/S0895479800376314. [8] D. Noutsos, P. Vassalos, Superlinear convergence for PCG using band plus algebra preconditioners for Toeplitz systems, Comput. Math. Appl. 56 (2008) 12551270, http://dx.doi.org/10.1016/j.camwa.2008.02.046. [9] D. Noutsos, S. Serra-Capizzano, P. Vassalos, Spectral equivalence and matrix algebra preconditioners for multilevel Toeplitz systems: a negative result, Contemp. Math. 323 (2003) 313–322, http://dx.doi.org/10.1090/conm/323/05711. [10] D. Noutsos, S. Serra-Capizzano, P. Vassalos, Matrix algebra preconditioners for multilevel Toeplitz systems do not insure optimal convergence rate, Theoret. Comput. Sci. 315 (2004) 557579, http://dx.doi.org/10.1016/j.tcs.2004.01.007. [11] D. Noutsos, S. Serra-Capizzano, P. Vassalos, Essential spectral equivalence via multiple step preconditioning and applications to ill conditioned Toeplitz matrices, Linear Algebra Appl. 491 (2016) 276–291, http://dx.doi.org/10.1016/j.laa.2015.08.021.

Please cite this article as: D. Noutsos and G. Tachyridis, Band Toeplitz preconditioners for non-symmetric real Toeplitz systems by preconditioned GMRES method, Journal of Computational and Applied Mathematics (2019), https://doi.org/10.1016/j.cam.2019.04.030.

18

D. Noutsos and G. Tachyridis / Journal of Computational and Applied Mathematics xxx (xxxx) xxx

[12] D. Bini, F. Di Benedetto, A new preconditioner for the parallel solution of positive definite toeplitz systems, in: Proceedings of the Second Annual ACM Symposium on Parallel Algorithms and Architectures, in: SPAA ’90, ACM, New York, NY, USA, 1990, pp. 220–223, http://dx.doi.org/10.1145/97444.97688. [13] F. Di Benedetto, S. Serra-Capizzano, A unifying approach to abstract matrix algebra preconditioning, Numer. Math. 82 (1999) 57–90, http: //dx.doi.org/10.1007/s002110050411. [14] X. Jin, Developments and Applications of Block Toeplitz Iterative Solvers, Science Press, Beijing, 2006. [15] S. Serra, Preconditioning strategies for asymptotically ill-conditioned block Toeplitz systems, BIT 34 (1994) 579–594, http://dx.doi.org/10.1137/ S0036142996311866. [16] M.K. Ng, Band preconditioners for block-Toeplitz-Toeplitz-block systems, Linear Algebra Appl. 259 (1997) 307–327, http://dx.doi.org/10.1016/ S0024-3795(96)00295-9. [17] D. Noutsos, S. Serra-Capizzano, P. Vassalos, A preconditioning proposal for ill-conditioned hermitian two-level Toeplitz systems, Numer. Linear Algebra Appl. 12 (2005) 231–239, http://dx.doi.org/10.1002/nla.398. [18] D. Noutsos, S. Serra-Capizzano, P. Vassalos, Block band Toeplitz preconditioners derived from generating function approximations: Analysis and applications, Numer. Math. 104 (2006) 339–376, http://dx.doi.org/10.1007/s00211-006-0020-7. [19] D. Noutsos, S. Serra-Capizzano, P. Vassalos, Two-level Toeplitz preconditioning: approximation results for matrices and functions, SIAM J. Sci. Comput. 28 (2006) 439–458, http://dx.doi.org/10.1137/050627046. [20] D. Noutsos, P. Vassalos, Band plus algebra preconditioners for two-level Toeplitz systems, BIT 51 (2011) 659–719, http://dx.doi.org/10.1007/ s10543-011-0314-8. [21] S. Serra-Capizzano, Matrix algebra preconditioners for multilevel Toeplitz matrices are not superlinear, Linear Algebra Appl. 343 (2002) 303319, http://dx.doi.org/10.1016/S0024-3795(01)00361-5. [22] S. Serra-Capizzano, E.E. Tyrtyshnikov, Any circulant-like preconditioner for multilevel matrices is not superlinear, SIAM J. Matrix Anal. Appl. 21 (2000) 431–439, http://dx.doi.org/10.1137/S0895479897331941. [23] S. Serra-Capizzano, E.E. Tyrtyshnikov, How to prove that a preconditioner cannot be superlinear, Math. Comp. 72 (2003) 13051316, http: //dx.doi.org/10.1090/S0025-5718-03-01506-0. [24] S.V. Parter, On the distribution of singular values of Toeplitz matrices, Linear Algebra Appl. 80 (1986) 115–130, http://dx.doi.org/10.1016/00243795(86)90280-6. [25] P. Tilli, Singular values and eigenvalues of non-Hermitian block Toeplitz matrices, Linear Algebra Appl. 272 (1998) 59–89, http://dx.doi.org/10. 1016/S0024-3795(97)00308-X. [26] P. Tilli, Some results on complex Toeplitz eigenvalues, J. Math. Anal. Appl. 239 (1999) 390–401, http://dx.doi.org/10.1006/jmaa.1999.6572. [27] E.E. Tyrtyshnikov, A unifying approach to some old and new theorems on distribution and clustering, Linear Algebra Appl. 232 (1996) 1–43, http://dx.doi.org/10.1016/0024-3795(94)00025-5. [28] E.E. Tyrtyshnikov, N.L. Zamarashkin, Thin structure of eigenvalue clusters for non-Hermitian Toeplitz matrices, Linear Algebra Appl. 292 (1999) 297–310, http://dx.doi.org/10.1016/S0024-3795(99)00044-0. [29] Y. Saad, Iterative Methods for Sparse Linear Systems, SIAM, Philadelphia, 2003, http://dx.doi.org/10.1137/1.9780898718003. [30] Y. Saad, M.H. Schultz, GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems, SIAM J. Sci. Stat. Comput. 7 (1986) 856–869, http://dx.doi.org/10.1137/0907058. [31] N.M. Nachtigal, S.C. Reddy, L.N. Trefethen, How fast are nonsymmetric matrix iterations? SIAM J. Matrix Anal. Appl. 13 (1992) 778795, http://dx.doi.org/10.1137/0613049. [32] E.E. Tyrtyshnikov, Influence of matrix operations on the distribution of eigenvalues and singular values of toeplitz matrices, Linear Algebra Appl. 207 (1994) 225–249, http://dx.doi.org/10.1016/0024-3795(94)90012-4. [33] H. Dai, Z. Geary, L.P. Kadanoff, Asymptotics of eigenvalues and eigenvectors of toeplitz matrices, J. Stat. Mech. Theory Exp. 2009 (05) (2009) P05012, http://dx.doi.org/10.1088/1742-5468/2009/05/P05012. [34] U. Grenander, G. Szegő, Toeplitz Forms and Their Applications, Chelsea, New York, 1984. [35] F. Avram, On bilinear forms in Gaussian random variables and toeplitz matrices, Probab. Theory Relat. Fields 79 (1988) 37–45, http: //dx.doi.org/10.1007/BF00319101. [36] S. Serra-Capizzano, D. Sesana, Tools for the eigenvalue distribution in a non-Hermitian setting, Linear Algebra Appl. 430 (2009) 423–437, http://dx.doi.org/10.1016/j.laa.2008.08.006. [37] M. Donatelli, M. Mazza, S. Serra-Capizzano, Spectral analysis and structure preserving preconditioners for fractional diffusion equations, J. Comput. Phys. 307 (2016) 262–279, http://dx.doi.org/10.1016/j.jcp.2015.11.061.

Please cite this article as: D. Noutsos and G. Tachyridis, Band Toeplitz preconditioners for non-symmetric real Toeplitz systems by preconditioned GMRES method, Journal of Computational and Applied Mathematics (2019), https://doi.org/10.1016/j.cam.2019.04.030.