Effective methods for solving banded Toeplitz systems

Dario Andrea Bini and Beatrice Meini Dipartimento di Matematica, Università di Pisa

            bini@dm.unipi.it, meini@dm.unipi.it

Abstract We propose new algorithms for solving n×n banded Toeplitz systems with bandwidth m. If the function associated with the Toeplitz matrix has no zero in the unit circle, then O(n log m + m log^{2} m log log ϵ^{−1}) arithmetic operations (ops) are sufficient to approximate the solution of the system up to within the error ϵ, otherwise the cost becomes O(n log m+m log^{2} m log \frac{n}{m}) ops. Here m = o(n), and n > log ϵ^{−1}. Some applications are presented. The methods can be applied to infinite and bi-infinite systems and to block matrices.

Key Words. Banded matrices, Toeplitz matrices, displacement rank, cyclic reduction, Graeffe’s method.

AMS Subject Classification. 65F05, 65F10, 15A23

1 Introduction

Let {a_{k}}_{k∈Z} be a sequence of real numbers such that a_{k} = 0 for k < m_{2}, and for k > m_{1}, a_{m_{1}}a_{m_{2}} ≠ 0, where Z is the ring of integers, and m_{1}, m_{2} ∈ Z, m_{2} < 0 < m_{1}. We may associate with {a_{k}}_{k∈Z} the function ϕ(z) = ∑_{i=m_{2}}^{m_{1}} z^{i}a_{i} and a matrix T = (t_{i,j}), having entries t_{i,j} = a_{i−j}. Let us use different notations for T according to its size. That is, we denote with T_{±∞} = (t_{i,j})_{i,j∈Z} the bi-infinite matrix, with T_{∞} = (t_{i,j})_{i,j≥0} the infinite matrix, and with T_{n} = (t_{i,j})_{i,j=1,...,n} the matrix having size n. We leave the notation T for a matrix having an unspecified size.

T is a banded matrix with bandwidth m = max{m_{1}, −m_{2}}, and has the Toeplitz structure, i.e., its entries are constant along the diagonals.

In this paper we propose a new approach to solving the (block) banded Toeplitz systems T x = b with bandwidth m, leading to an algorithm that can be used either as a direct algorithm or as an iterative one. For n × n systems the direct method, which can be applied with almost no restriction, has a sequential cost of O(n log m + t(m) log \frac{n}{m}) ops, where t(m) is the arithmetic cost of solving an m × m Toeplitz-like system, i.e., t(m) ∈ {m^{2}, m log^{2} m, m log m} according to the algorithm adopted (compare [23], [1], [15]).

Under suitable conditions generally satisfied in most part of the applications, the direct algorithm turns into an iterative one, and its convergence speed is quadratic. In this way the number of ops needed to approximate n components of the solution up to within an O(ϵ) error is O(n log m + t(m) log log ϵ^{−1}).

The cost O(t(m) log log ϵ^{−1}) is required for arriving at an approximate LU factorization of the (infinite) matrix obtained by means of a suitable permutation of rows and columns of T. The O(n log m) term is the cost of the back substitution stage of the algorithm. This nice convergence property allows us to approximate n consecutive components of the solution of infinite and bi-infinite systems with the same computational cost as well. In the PRAM model of parallel computation, the algorithm costs O(t_{s}(m) log \frac{n}{m}) steps with max\{t_{p}(m), n\} processors, where t_{s}(m), t_{p}(m) are the number of steps and of processors, respectively, needed to solve an m × m Toeplitz-like system.

1.1 The problem and its relevance.

The problem of solving banded Toeplitz systems as T x = b, plays an important role in many applications. Example of these systems are frequently encountered in the numerical solution of ordinary differential equations with boundary values, where the banded Toeplitz matrix provides a finite difference discretization of the differential operator [31], and in the numerical solution of many Markov chains coming from the modeling of queueing problems, where the involved matrix defines the probabilities of transition from one state to the other [26], [25], [18]. Another indirect and more recent application is related to the solution of general Toeplitz systems, like the ones encountered in the problems in signal processing [16], where the banded Toeplitz matrix is used as a preconditioner for speeding up the convergence of preconditioned conjugate gradient techniques [17], [30], [29].

In all the above problems the dimension of the matrices involved is typically huge or, sometimes, may be even infinite [26]. Moreover, in multi dimensional problems the banded Toeplitz structure is expressed in terms of blocks, i.e., the scalar entries a_{j−i} are replaced by square matrices.

This makes it particularly important to devise efficient and reliable algorithms for the solution of (block) banded (block) Toeplitz systems that maintain their efficiency in the case of matrices having a very large, or even infinite, size.

A certain attention to this problem has been devoted in the literature specially concerning the analysis of the complexity in parallel models of computation, like the PRAM model. In the PRAM model it is assumed that a certain number of processors work concurrently, more specifically, at each time unit (parallel step) each processor performs a single arithmetic operation. In this model the cost of a computation is expressed in terms of the number of parallel steps and of the number of processors needed to carry out the computation.

Different direct parallel algorithms have been developed [20], [4], [6], [5]. They are based mainly on two different techniques: embedding the matrix into a slightly larger matrix that can be diagonalized by means of trigonometric transforms (FFT); correcting the matrix by means of a low rank additive correction in order to find a new matrix that can be diagonalized by a trigonometric transform and applying the Morrison-Sherman-Woodbury formula.

The algorithms obtained in this way have a parallel cost of O(log n + log^{2} m) steps with nm^{2} processors, and their sequential cost of O(n log n + m^{2}n) arithmetic operations (ops) is not inferior to the cost of Gaussian elimination, that is, O(m^{2}n) ops. Moreover most of these algorithms have problems of numerical stability.

Recently, a direct and fully parallelizable algorithm for inverting block Toeplitz matrices in block Hessenberg form, based on a doubling technique and on the use of the Morrison-Sherman-Woodbury formula, has been proposed in [11]. The algorithm, applied to solve the system T_{n}x = b, has a sequential cost of O(m^{2} log \frac{n}{m} + n log n) ops. However, this algorithm does not fully exploit the band structure of the system and requires a large storage.

Iterative algorithms, based on the preconditioned conjugate gradient (pcg) techniques, are fully parallelizable and have a sequential cost per iteration of O(n log n). There exists a wide literature concerning these techniques, see for instance [16], [15], [29], [32]. The number of pcg iterations needed in the worst case depends on the bandwidth m (compare [30] for the dependency of the number of iterations on the specific features of the function associated with the matrix). Due to the Axelsson-Lindskog theorem [2], the number of iterations sufficient to reach a super-linear convergence is proportional to m if a circulant [14], [32] or τ-class [6], [7], [28] preconditioner is used. Moreover, O(mn log n) ops are needed for the presently known preconditioners [30].

1.2 New results.

The new approach is based on the following idea. The matrix T is partitioned into m×m blocks, thus obtaining a block Toeplitz block tridiagonal matrix with Toeplitz blocks. The technique of cyclic reduction (see [19]) is applied to the reblocked matrix T, that is, an odd-even permutation of the block rows and the block columns of T is followed by a step of block Gaussian elimination.

We prove that the block tridiagonal block Toeplitz structure is kept almost everywhere by the cyclic reduction step. In other words, the Schur complement that we obtain after applying the single step of block Gaussian elimination to the matrix with permuted block rows and block columns is still a block tridiagonal block Toeplitz matrix except for the block in position (1,1).

On the other hand a sort of deterioration occurs in the inner Toeplitz structure of the blocks of T. In fact, the blocks of the Schur complement obtained just at the first step of cyclic reduction are not Toeplitz any more. In principle, this fact would imply that each block must be stored in m^{2} memory locations and that managing with blocks would cost O(m^{3}) ops. In this way cyclic reduction would require O(m^{2}n + m^{3} log \frac{n}{m}) ops in order to solve an n × n banded Toeplitz system with bandwidth m.

Even though the Toeplitz structure of the blocks is lost, a more general structure is maintained by the blocks of the Schur complement at each level of cyclic reduction. More precisely, the main result that we prove in section 3 is that each block defining the Schur complement is a Toeplitz-like matrix of displacement rank at most four. Each block can be represented as the sum of at most four products of pairs of lower and upper triangular Toeplitz matrices. Therefore it can be stored in O(m) memory locations by representing the block by means of few vectors. Furthermore, performing multiplications between the blocks costs O(m log m) ops., and inverting a block costs t(m) ∈ {m^{2}, m log^{2} m, m log m} ops, t(m) ∈ {m^{2}, m log^{2} m, m log m}.

Moreover we prove explicit relations which hold among the vectors defining the blocks of the Schur complement at two consecutive steps. These relations allow us to implement the cyclic reduction technique in O(n log m + t(m) log \frac{n}{m}) ops.

The algorithm can be easily extended to banded block Toeplitz systems. For blocks of size h the cost is increased by the factor h^{3}.

The cyclic reduction algorithm obtained in this way has a further nice feature. If the polynomial a(z) = z^{−m_{2}}ϕ(z) = ∑_{i=m_{2}}^{m_{1}} a_{i}z^{i−m_{2}} has m zeros inside (outside) the unit circle and if the blocks have bounded norm, we may prove that the lower (upper) diagonal blocks of the Schur complement obtained at the j-th step of cyclic reduction, tend to zero as γ^{2^{j}}, where 0 < γ < 1 (double exponential convergence).

For stochastic matrices, like the ones modeling Markov chains that arise in QBD problems [26], the above conditions are satisfied [9], [10]. In the case where T is symmetric and the polynomial a(z) has no zeros of modulus 1, then the above property is satisfied simultaneously inside and outside the unit circle, and the boundness of the blocks follows from the Cauchy interlace theorem. In this case we may prove that the constant γ, that rules the double exponential convergence, is the modulus of the zero of a(z) in the open unit disc, closest to the unit circle.

This property generalizes and extends the results of [13] proved for tridiagonal Toeplitz matrices, where m = 1. The conditions under which the convergence occurs are much weaker than the conditions given in [22], [33].

For matrices associated with polynomials having some zeros of modulus 1, we provide a simple trick that enables us to move the zeros of modulus 1 slightly in/out the unit disc. Then, constructing the normal equations leads us to a system having a banded matrix with a nearly Toeplitz structure, that satisfies the convergence theorem. This fact allows us to improve the robustness of our algorithm.

The algorithm proposed in this paper is numerically stable if applied to positive definite symmetric matrices. In fact, all the m × m linear systems involved at each step of cyclic reduction have condition numbers not greater than the condition number of the original matrix T.

Another interesting feature of the algorithm is its close relationship with polynomial computations. In fact, we show that applying the cyclic reduction algorithm to a banded Toeplitz matrix corresponds to applying the Graeffe iteration [27] to the polynomial α(z) having as zeros the m-th powers of the zeros of a(z). Moreover, this similarity between polynomial and Toeplitz computations enables us to extend the Graeffe iteration to matrix polynomials.

We have implemented in Fortran 90 different versions of our algorithm which include the cases of matrices with scalar/block entries, symmetric/unsymmetric matrices, and stochastic matrices.

The algorithms have been compared with the preconditioned conjugate gradient method with circulant preconditioner. From the many numerical experiments performed over an alpha workstation, our algorithms have proved to be robust, numerically stable and fast. For matrices of size n = 524288 and bandwidth m = 32 our programs have been 16 times faster than the preconditioned conjugate gradient method. For larger dimensions the conjugate gradient method could not be applied for lack of memory.

The programs can be taken by anonymous ftp at the following URL address: morse.dm.unipi.it.

2 Preliminaries

Without loss of generality let us assume that −m_{2} ≥ m_{1}, moreover, let n be an integer multiple of m, i.e., n = mq, for an integer q. We refer to ϕ(z) as the generating function of the matrix T and we call T_{n} an n-section of T_{∞}.

The matrix T can be partitioned into m × m blocks yielding a block tridiagonal Toeplitz matrix. In particular T_{n} has the following structure

T_{n} = [ \begin{matrix}A_{0} & A_{−1} & & & ◯ \\ A_{1} & A_{0} & A_{−1} \\ & ... & ... & ... \\ & & A_{1} & A_{0} & A_{−1} \\ ◯ & & & A_{1} & A_{0} \end{matrix} ]

where

\begin{matrix}A_{1} = [ \begin{matrix}a_{m} & a_{m−1} & … & a_{1} \\ 0 & a_{m} & ... & ⋮ \\ ⋮ & ... & ... & a_{m−1} \\ 0 & … & 0 & a_{m} \end{matrix} ] , \\ A_{0} = [ \begin{matrix}a_{0} & a_{−1} & … & a_{−m+1} \\ a_{1} & a_{0} & ... & ⋮ \\ ⋮ & ... & ... & a_{−1} \\ a_{m−1} & … & a_{1} & a_{0} \end{matrix} ] \\ A_{−1} = [ \begin{matrix}a_{−m} & 0 & … & 0 \\ a_{−m+1} & a_{−m} & ... & ⋮ \\ ⋮ & ... & ... & 0 \\ a_{−1} & … & a_{−m+1} & a_{−m} \end{matrix} ] \end{matrix}  (1)

and similarly T_{∞} and T_{±∞}.

In the light of the above block structure we may associate with T the matrix polynomial A(z) = A_{−1} + zA_{0} + z^{2}A_{1}.

In this section we first recall the cyclic reduction algorithm, originally introduced for the numerical solution of the Poisson equation ([19]), used in [8] for computing the probability invariant vector of stochastic matrices, and here adjusted in order to solve the block tridiagonal block Toeplitz system T x = b. Then we introduce the concept of displacement rank and we recall some properties that are needed for the design and analysis of our algorithm.

2.1 Cyclic reduction.

Consider the system

T x = b  (2)

and partition the vectors x and b into blocks x_{k}, b_{k} of dimension m, respectively, that is, x = (x_{k}), b = (b_{k}). By performing an odd-even permutation of the block rows and block columns in (2) we find that

[ \begin{matrix}D_{1}^{(0)} & U^{(0)} \\ L^{(0)} & D_{2}^{(0)} \end{matrix} ] [ \begin{matrix}x & {(0) \atop +} \\ x & {(0) \atop −} \end{matrix} ] = [ \begin{matrix}b & {(0) \atop +} \\ b & {(0) \atop −} \end{matrix} ]  (3)

where x_{+}^{(0)} = (x_{2k}), x_{−}^{(0)} = (x_{2k−1}), b_{+}^{(0)} = (b_{2k}), b_{−}^{(0)} = (b_{2k−1});

\begin{matrix}D_{1}^{(0)} = D_{2}^{(0)} = [ \begin{matrix}A_{0} & & ◯ \\ & A_{0} \\ ◯ & & ... \end{matrix} ] , L^{(0)} = [ \begin{matrix}A_{−1} & & & ◯ \\ A_{1} & A_{−1} \\ & ... & ... \\ ◯ & & ... & ... \end{matrix} ] , \\ U^{(0)} = [ \begin{matrix}A_{1} & A_{−1} & & ◯ \\ & A_{1} & A_{−1} \\ & & ... & ... \\ ◯ & & & ... \end{matrix} ] . \end{matrix}

Applying one step of block-Gaussian elimination to the 2 × 2 block system (3) yields the equivalent system

\{ \begin{matrix}(D_{2}^{(0)} − L^{(0)}D_{1}^{(0)^{−1}}U^{(0)})x_{−}^{(0)} = b_{−}^{(0)} − L^{(0)}D_{1}^{(0)^{−1}}b_{+}^{(0)} \\ x_{+}^{(0)} = D_{1}^{(0)}^{−1}(b_{+}^{(0)} − U^{(0)}x_{−}^{(0)}). \end{matrix}  (4)

Denoting

T^{(1)} = D_{2}^{(0)} − L^{(0)}D_{1}^{(0)^{−1}}U^{(0)}  (5)

the Schur complement of D_{2}^{(0)}, x^{(1)} = x_{−}^{(0)}, b^{(1)} = b_{−}^{(0)} − L^{(0)}D_{1}^{(0)^{−1}}b_{+}^{(0)}, the above system is ultimately reduced to solving

T^{(1)}x^{(1)} = b^{(1)}.  (6)

It is interesting to observe that T^{(1)} is a block tridiagonal matrix which, except for the north west corner block, has the block Toeplitz structure, i.e.,

T^{(1)} = [ \begin{matrix}Â^{(1)} & A_{−1}^{(1)} & & ◯ \\ A_{1}^{(1)} & A_{0}^{(1)} & A_{−1}^{(1)} \\ & A_{1}^{(1)} & A_{0}^{(1)} & ... \\ ◯ & & ... & ... \end{matrix} ] .

In other words, T^{(1)} is uniquely determined by the blocks A_{−1}^{(1)}, A_{0}^{(1)}, A_{1}^{(1)}, Â^{(1)}. Once x^{(1)} has been computed or approximated, the solution of the original system (2) can be recovered by means of back substitution through (4). In order to compute the solution x^{(1)} we may cyclically apply the same reduction (cyclic reduction) to the system (6). In this way we obtain a sequence of systems {T^{(j)}x^{(j)} = b^{(j)}}, where T^{(0)} = T, and T^{(j)} is a block tridiagonal matrix having, except for the north west corner, the block Toeplitz structure. Each matrix T^{(j)} is uniquely determined by the blocks Â^{(j)}, A_{−1}^{(j)}, A_{0}^{(j)}, A_{1}^{(j)}.

For a finite banded Toeplitz system T_{n}x = b, where n = mq, q = 2^{p}, p positive integer, the cyclic reduction generates a sequence of systems T_{n_{j}}^{(j)}x^{(j)} = b^{(j)}, j = 1, 2, …, p, of dimension n_{j} = m2^{p−j}.

For a bi-infinite block tridiagonal block Toeplitz system T_{±∞}x = b, where x = (x_{i})_{i∈Z}, b = (b_{i})_{i∈Z}, the cyclic reduction generates a sequence of bi-infinite block tridiagonal block Toeplitz systems T_{±∞}^{(j)}x^{(j)} = b^{(j)} where T_{±∞}^{(j)} is the bi-infinite block tridiagonal block Toeplitz matrix defined by the m × m blocks A_{i}^{(j)}, i = −1, 0, 1.

For finite matrices of size n = m2^{p}, the cyclic reduction process is carried out in p steps. The solution x^{(p)} is computed by solving the m×m system T^{(p)}x^{(p)} = b^{(p)} and the whole solution of (2) is computed by means of back substitution, recursively applying (4) extended to the generic j-th step.

For infinite or bi-infinite matrices, under suitable conditions (see section 4), we may prove that the sequence of systems T^{(j)}x^{(j)} = b^{(j)} converges to a system having a block (bi-)diagonal structure. In this way, by choosing j large enough, we may approximate to any precision an arbitrary finite number of components of the solution x^{(j)}. The back substitution stage completes the recovering of the solution of (2).

The block entries A_{i}^{(j+1)}, i = −1, 0, 1, Â^{(j+1)}, obtained at step j + 1 that define the matrix T^{(j+1)}, are related to the block entries A_{i}^{(j)}, i = −1, 0, 1, Â^{(j)}, obtained at step j, by means of simple relations. For this purpose let us associate with T^{(j)} the following matrix polynomial

A^{(j)}(z) = A_{−1}^{(j)} + zA_{0}^{(j)} + z^{2}A_{1}^{(j)}  (7)

and express in compact form the block odd-even permutation and the Gaussian elimination of one step of cyclic reduction in terms of A^{(j)}(z).

In fact from (5) we may easily deduce that

\{ \begin{matrix}A^{(j+1)}(z) = zA_{0}^{(j)} − (A_{−1}^{(j)} + zA_{1}^{(j)})A_{0}^{(j)^{−1}}(A_{−1}^{(j)} + zA_{1}^{(j)}) \\ Â^{(j+1)} = Â^{(j)} − A_{−1}^{(j)}A_{0}^{(j)^{−1}}A_{1}^{(j)}. \end{matrix}  (8)

The above relations yield the following matrix equations

\{ \begin{matrix}A_{−1}^{(j+1)} = −A_{−1}^{(j)}A_{0}^{(j)^{−1}}A_{−1}^{(j)} \\ A_{0}^{(j+1)} = A_{0}^{(j)} − A_{1}^{(j)}A_{0}^{(j)^{−1}}A_{−1}^{(j)} − A_{−1}^{(j)}A_{0}^{(j)^{−1}}A_{1}^{(j)} \\ A_{1}^{(j+1)} = −A_{1}^{(j)}A_{0}^{(j)^{−1}}A_{1}^{(j)} \\ Â^{(j+1)} = Â^{(j)} − A_{−1}^{(j)}A_{0}^{(j)^{−1}}A_{1}^{(j)} \end{matrix}  (9)

which allow us to compute the sequence of blocks defining the matrices T^{(j)}, j = 1, 2, . . ..

A direct inspection shows that the Toeplitz structure of the blocks A_{i}^{(0)}, i = −1, 0, 1, is not kept any more by the blocks A_{i}^{(j)}, i = −1, 0, 1, and Â^{(j)}. For this property one might deduce that the initial structure cannot be exploited in order to increase the efficiency of the cyclic reduction algorithm.

However, we may prove that a more general structure, i.e., the displacement structure is maintained by the blocks A_{i}^{(j)}, i = −1, 0, 1, and Â^{(j)}. This fact allows us to devise a suitable FFT-based implementation of the cyclic reduction algorithm thus reaching a dramatic reduction of the complexity bound.

For this purpose we need to rewrite the functional relation (8) for the matrix power series A^{(j)}(z) in the following equivalent form:

A^{(j+1)}(z^{2}) = (\frac{A^{(j)}(z)^{−1} − A^{(j)}(−z)^{−1}}{2z} )^{−1} = −A^{(j)}(z)A_{0}^{(j)^{−1}}A^{(j)}(−z). (10)

2.2 Displacement rank.

The concept of displacement rank, introduced by Kailath et al. in [24], is fundamental in devising and analyzing algorithms related to Toeplitz matrices. Several displacement operators have been introduced and analyzed in the literature; here we adopt an operator that seems to be particularly suitable for dealing with our problems.

Define the m × m down-shift matrix

Z = [ \begin{matrix}0 \\ 1 & 0 \\ & ... & ... \\ & & 1 & 0 \end{matrix} ] ,

and consider the following displacement operator

Δ(H) = ZH − HZ

defined for any m × m matrix H. Moreover, denote by L(w) the m × m lower triangular Toeplitz matrix defined by its first column w.

Let us also introduce the following notations: the vectors e_{1} = [1, 0, …, 0]^{T}, e_{m} = [0, …, 0, 1]^{T} denote the first and the last column of the identity matrix of order m, respectively.

Observe that for a general Toeplitz matrix A the displacement Δ(A) is zero except for the entries in the first row and in the last column, whence the rank of Δ(A) is at most 2.

We say that the matrix H has displacement rank r if r is the rank of Δ(H), i.e., if there exist vectors u_{i}, v_{i}, i = 1, …, r, satisfying the equation Δ(H) = ∑_{i=1}^{r} u_{i}v_{i}^{T}.

According to this definition an m × m Toeplitz matrix A has displacement rank at most 2, more precisely it holds

Δ(A) = −e_{1}e_{m}^{T}AZ + ZAe_{m}e_{m}^{T}.

A nice property of Δ(H) is that the displacement of the inverse matrix can be explicitly related to the displacement of the matrix itself, in fact, if H is nonsingular, then

Δ(H^{−1}) = −H^{−1}Δ(H)H^{−1}.  (11)

The following result has been proved in [3], [12].

Theorem 1 Let K be an m × m matrix such that Δ(K) = ∑_{i=1}^{r} u_{i}v_{i}^{T} where u_{i} and v_{i} are m-dimensional vectors. Then we have

K = L(Ke_{1}) − ∑_{i=1}^{r} L(u_{i})L^{T} (Zv_{i}).

The above result allows one to represent any matrix K as a sum of products of lower and upper triangular Toeplitz matrices defined by the first column of K and by the vectors u_{i}, v_{i} associated with the displacement of K.

If the matrix K is nonsingular then the above representation theorem can be applied to K^{−1} in the light of (11) yielding the following result

K^{−1} = L(K^{−1}e_{1}) + ∑_{i=1}^{r} L(K^{−1}u_{i})L^{T} (ZK^{−T} v_{i}).

We recall that multiplying an m×m Toeplitz matrix and a vector costs O(m log m) ops if the computation is performed by means of FFT [12]. In the PRAM model, the cost is O(log m) steps with m processors.

3 Structural properties of cyclic reduction.

In this section we prove that the matrix polynomials A^{(j)}(z), j = 0, 1, …, defined in (7), (10), have displacement rank at most 2, and provide explicit formulae relating the displacement of A^{(j)}(z) to the displacement of A^{(j+1)}(z). As a consequence we find that A_{−1}^{(j)} and A_{1}^{(j)} have displacement rank at most 2, A_{0}^{(j)} has displacement rank at most 4. Moreover, we show that the blocks Â^{(j)} have displacement rank at most 3.

These properties allow us to devise an algorithm, based on FFT, that performs a single cyclic reduction step in O(m log^{2} m) arithmetic operations if super-fast algorithms for the inversion of m × m Toeplitz-like matrices are used.

We need the following result.

Theorem 2 For the matrix polynomial A^{(j)}(z) generated by the cyclic reduction (10) we have

Δ(A^{(j)}(z)) = a^{(j)}(z)u^{(j)T} (z) − v^{(j)}(z)c^{(j)T} (z), j ≥ 0,  (12)

where the vectors a^{(j)}(z), c^{(j)}(z), u^{(j)}(z), v^{(j)}(z) have entries which are polynomials in z, and are defined by the following equations:

\begin{matrix}a^{(j)}(z) = A^{(j)}(z)e_{1}, c^{(j)T} (z) = e_{m}^{T}A^{(j)}(z), j ≥ 0 \\ u^{(j)}(z^{2})^{T} = (−u^{(j−1)}(z)^{T} (A_{0}^{(j−1)})^{−1}A^{(j−1)}(−z)+ \\ u^{(j−1)}(−z)^{T} (A_{0}^{(j−1)})^{−1}A^{(j−1)}(z))/(2z), j ≥ 1; \\ v^{(j)}(z^{2}) = (−A^{(j−1)}(−z)(A_{0}^{(j−1)})^{−1}v^{(j−1)}(z)+ \\ A^{(j−1)}(z)(A_{0}^{(j−1)})^{−1}v^{(j−1)}(−z))/(2z), j ≥ 1; \\ u^{(0)}(z)^{T} = ze_{m}^{T}, v^{(0)}(z) = ze_{1}. \end{matrix} (13)

Proof. We prove relation (12) by induction on j. For j = 0 a simple calculation shows that

\begin{matrix}Δ(A^{(0)}(z)) = z(A^{(0)}(z)e_{1}e_{m}^{T} − e_{1}e_{m}^{T}A^{(0)}(z)) = \\ A^{(0)}(z)e_{1}u^{(0)T} (z) − v^{(0)}(z)e_{m}^{T}A^{(0)}(z). \end{matrix}

Suppose that relation (12) holds for a fixed j ≥ 0; we show that it holds for j + 1. For simplicity let us denote B = A^{(j)}(z)^{−1} − A^{(j)}(−z)^{−1}. For the properties of the displacement operator (11) and for (10) we have:

\begin{matrix}Δ(A^{(j+1)}(z^{2})) = 2zΔ(B^{−1}) = −2zB^{−1}Δ(B)B^{−1} = \\ −2zB^{−1}(−A^{(j)}(z)^{−1}Δ(A^{(j)}(z))A^{(j)}(z)^{−1} + A^{(j)}(−z)^{−1}Δ(A^{(j)}(−z))A^{(j)}(−z)^{−1})B^{−1}; \end{matrix}

by substituting relation (12) in the above equation we obtain

Δ(A^{(j+1)}(z^{2})) = 2zB^{−1}e_{1}u^{(j)}(z)^{T} A^{(j)}(z)^{−1}B^{−1} − 2zB^{−1}A^{(j)}(z)^{−1}v^{(j)}(z)e_{m}^{T}B^{−1}− 2zB^{−1}e_{1}u^{(j)}(−z)^{T} A^{(j)}(−z)^{−1}B^{−1} + 2zB^{−1}A^{(j)}(−z)^{−1}v^{(j)}(−z)e_{m}^{T}B^{−1} = A^{(j+1)}(z^{2})e_{1}u^{(j+1)}(z^{2})^{T} − v^{(j+1)}(z^{2})e_{m}^{T}A^{(j+1)}(z^{2}).

□

By induction on j, from (13), it readily follows that u^{(j)}(z) and v^{(j)}(z) have components which are polynomials of degree at most 1. Therefore, since a^{(j)}(z) and c^{(j)}(z) have degree at most 2, we may denote

u^{(j)}(z) = u_{−1}^{(j)} + zu_{0}^{(j)},
v^{(j)}(z) = v_{−1}^{(j)} + zv_{0}^{(j)},
a^{(j)}(z) = a_{−1}^{(j)} + za_{0}^{(j)} + z^{2}a_{1}^{(j)},
c^{(j)}(z) = c_{−1}^{(j)} + zc_{0}^{(j)} + z^{2}c_{1}^{(j)}.

Moreover, by comparing the term of degree 3 in (12) we may deduce that, for any j ≥ 0, it holds

a_{1}^{(j)}u_{0}^{(j)T} = v_{0}^{(j)}c_{1}^{(j)T} .  (14)

In this way, from the above properties and from Theorem 2 we deduce explicit relations for the vectors u_{i}^{(j)}, v_{i}^{(j)}, i = −1, 0, a_{i}^{(j)}, c_{i}^{(j)}, i = −1, 0, 1, j ≥ 0, defining Δ(A^{(j)}(z)):

\begin{matrix}u_{−1}^{(j+1)T} = u_{−1}^{(j)T} − u_{0}^{(j)T} A_{0}^{(j)^{−1}}A_{−1}^{(j)} \\ u_{0}^{(j+1)T} = −u_{0}^{(j)T} A_{0}^{(j)^{−1}}A_{1}^{(j)} \\ v_{−1}^{(j+1)} = v_{−1}^{(j)} − A_{−1}^{(j)}A_{0}^{(j)^{−1}}v_{0}^{(j)} \\ v_{0}^{(j+1)} = −A_{1}^{(j)}A_{0}^{(j)^{−1}}v_{0}^{(j)} \\ a_{−1}^{(j+1)} = −A_{−1}^{(j)}A_{0}^{(j)^{−1}}a_{−1}^{(j)} \\ a_{0}^{(j+1)} = a_{0}^{(j)} − A_{−1}^{(j)}A_{0}^{(j)^{−1}}a_{1}^{(j)} − A_{1}^{(j)}A_{0}^{(j)^{−1}}a_{−1}^{(j)} \\ a_{1}^{(j+1)} = −A_{1}^{(j)}A_{0}^{(j)^{−1}}a_{1}^{(j)} \\ c_{−1}^{(j+1)T} = −c_{−1}^{(j)^{T}} A_{0}^{(j)^{−1}}A_{−1}^{(j)} \\ c_{0}^{(j+1)T} = c_{0}^{(j)T} − c_{−1}^{(j)T} A_{0}^{(j)^{−1}}A_{1}^{(j)} − c_{1}^{(j)T} A_{0}^{(j)^{−1}}A_{−1}^{(j)} \\ c_{1}^{(j+1)T} = −c_{1}^{(j)T} A_{0}^{(j)^{−1}}A_{1}^{(j)} \end{matrix}  (15)

where u_{−1}^{(0)} = v_{−1}^{(0)} = 0, u_{0}^{(0)} = e_{m}, v_{0}^{(0)} = e_{1}.

The above relations, together with Theorem 2, yield the following result:

Theorem 3 For the matrices A_{−1}^{(j)}, A_{0}^{(j)}, A_{1}^{(j)}, j ≥ 0, generated by the cyclic reduction (8) we have

\begin{matrix}Δ(A_{−1}^{(j)}) = a_{−1}^{(j)}u_{−1}^{(j)T} − v_{−1}^{(j)}c_{−1}^{(j)T} \\ Δ(A_{0}^{(j)}) = a_{−1}^{(j)}u_{0}^{(j)T} + a_{0}^{(j)}u_{−1}^{(j)T} − v_{−1}^{(j)}c_{0}^{(j)T} − v_{0}^{(j)}c_{−1}^{(j)T} \\ Δ(A_{1}^{(j)}) = r^{(j)}u_{0}^{(j)T} − v_{0}^{(j)}ĉ^{(j)T} , \end{matrix}  (16)

where

\begin{matrix}r^{(0)} = a_{0}^{(0)}, \\ r^{(j+1)} = r^{(j)} − A_{1}^{(j)}A_{0}^{(j)^{−1}}a_{−1}^{(j)}, j ≥ 0, \end{matrix}  (17)

and ĉ^{(j)T} = e_{m}^{T}Â^{(j)}.

Proof. By comparing the terms in z of degree 0 and 1 in relation (12) we obtain the expression for Δ(A_{−1}^{(j)}) and Δ(A_{0}^{(j)}), respectively. Similarly, by comparing the terms of degree 2 we have

Δ(A_{1}^{(j)}) = a_{0}^{(j)}u_{0}^{(j)T} + a_{1}^{(j)}u_{−1}^{(j)T} − v_{−1}^{(j)}c_{1}^{(j)T} − v_{0}^{(j)}c_{0}^{(j)T} .  (18)

For j = 0 relation (16) immediately follows. For j ≥ 1, from an inductive application of (15) we may deduce that u_{−1}^{(j)T} = − ∑_{i=0}^{j−1} u_{0}^{(i)T} A_{0}^{(i)^{−1}}A_{−1}^{(i)} and v_{−1}^{(j)} =

− ∑_{i=0}^{j−1} A_{−1}^{(i)}A_{0}^{(i)^{−1}}v_{0}^{(i)}, whence

a_{1}^{(j)}u_{−1}^{(j)T} = − ∑_{i=0}^{j−1} a_{1}^{(j)}u_{0}^{(i)T} A_{0}^{(i)^{−1}}A_{−1}^{(i)}.  (19)

On the other hand, from (9), it can be easily proved by induction on j that, for any j ≥ 1 and 0 ≤ i < j, it holds that

A_{1}^{(j)} = (−1)^{j−i} (∏_{h=i}^{j−1} A_{1}^{(j−h+i−1)}A_{0}^{(j−h+i−1)^{−1}}) A_{1}^{(i)},  (20)

where ∏_{h=h_{1}}^{h_{2}} X_{σ(h)} = X_{σ(h_{1})} ⋯ X_{σ(h_{2})}. Whence, from (19), (20) and (14), we have

\begin{matrix}a_{1}^{(j)}u_{−1}^{(j)T} = − ∑_{i=0}^{j−1}(−1)^{j−i} (∏_{h=i}^{j−1} A_{1}^{(j−h+i−1)}A_{0}^{(j−h+i−1)^{−1}}) a_{1}^{(i)}u_{0}^{(i)T} A_{0}^{(i)^{−1}}A_{−1}^{(i)} = \\ − ∑_{i=0}^{j−1}(−1)^{j−i} (∏_{h=i}^{j−1} A_{1}^{(j−h+i−1)}A_{0}^{(j−h+i−1)^{−1}}) v_{0}^{(i)}c_{1}^{(i)T} A_{0}^{(i)^{−1}}A_{−1}^{(i)}. \end{matrix}

(21)

By applying the recursive relation for v_{0}^{(j)}, given in formula (15), it may be proved by induction on j that, for any j ≥ 1, 0 ≤ i < j, it holds that

v_{0}^{(j)} = (−1)^{j−i} (∏_{h=i}^{j−1} A_{1}^{(j−h+i−1)}A_{0}^{(j−h+i−1)^{−1}}) v_{0}^{(i)}  (22)

Whence, from (21) and (22), it follows that

a_{1}^{(j)}u_{−1}^{(j)T} = − ∑_{i=0}^{j−1} v_{0}^{(j)}c_{1}^{(i)T} A_{0}^{(i)^{−1}}A_{−1}^{(i)},

and, in equation (18), we obtain

a_{1}^{(j)}u_{−1}^{(j)T} − v_{0}^{(j)}c_{0}^{(j)T} = −v_{0}^{(j)} (c_{0}^{(j)T} + ∑_{i=0}^{j−1} c_{1}^{(i)T} A_{0}^{(i)^{−1}}A_{−1}^{(i)}) = −v_{0}^{(j)}ĉ^{(j)T} .

Similarly, it may be proved that

a_{0}^{(j)}u_{0}^{(j)T} − v_{−1}^{(j)}c_{1}^{(j)T} = (a_{0}^{(j)} + ∑_{i=0}^{j−1} A_{−1}^{(i)}A_{0}^{(i)^{−1}}a_{1}^{(i)}) u_{0}^{(j)T} = r^{(j)}u_{0}^{(j)T}

□

The above theorem, in the light of the results of Section 2.2, yields the following representations for the matrices A_{−1}^{(j)}, A_{0}^{(j)}, A_{1}^{(j)}:

Theorem 4 At each step j of cyclic reduction, the matrices A_{−1}^{(j)}, A_{0}^{(j)}, A_{1}^{(j)}, A_{0}^{(j)^{−1}} can be rewritten as

A_{−1}^{(j)} = L(a_{−1}^{(j)})(I − L^{T} (Zu_{−1}^{(j)}) + L(v_{−1}^{(j)})L^{T} (Zc_{−1}^{(j)}))

A_{0}^{(j)} = L(a_{0}^{(j)})(I − L^{T} (Zu_{−1}^{(j)})) − L(a_{−1}^{(j)})L^{T} (Zu_{0}^{(j)}) + L(v_{−1}^{(j)})L^{T} (Zc_{0}^{(j)})+

L(v_{0}^{(j)})L^{T} (Zc_{−1}^{(j)})

A_{1}^{(j)} = L(a_{1}^{(j)}) − L(r^{(j)})L^{T} (Zu_{0}^{(j)}) + L(v_{0}^{(j)})L^{T} (Zĉ^{(j)})

A_{0}^{(j)^{−1}} = L(A_{0}^{(j)^{−1}}e_{1}) + L^{T} (ZA_{0}^{(j)^{−T}} u_{−1}^{(j)}) + L(A_{0}^{(j)^{−1}}a_{−1}^{(j)})L^{T} (ZA_{0}^{(j)^{−T}} u_{0}^{(j)})− L(A_{0}^{(j)^{−1}}v_{0}^{(j)})L^{T} (ZA_{0}^{(j)^{−T}} c_{−1}^{(j)}).

Proof. The proof readily follows from Theorem 3 and from Theorem 1. □ Similar properties hold for the matrices Â^{(j)}:

Theorem 5 For the matrices Â^{(j)}, j ≥ 0, generated by the cyclic reduction (8) we have

Δ(Â^{(j)}) = −e_{1}c_{−1}^{(0)T} + a_{−1}^{(j)}u_{0}^{(j)T} − v_{−1}^{(j)}ĉ^{(j)T} ,  (23)

whence the matrices Â^{(j)} can be rewritten as

Â^{(j)} = L(â^{(j)}) + L^{T} (Zc_{−1}^{(0)}) − L(a_{−1}^{(j)})L^{T} (Zu_{0}^{(j)}) + L(v_{−1}^{(j)})L^{T} (Zĉ^{(j)}),

where â^{(j)} = Â^{(j)}e_{1}, ĉ^{(j)T} = e_{m}^{T}Â^{(j)}.

Proof. We prove (23) by induction on j. For j = 0 the equation is obviously verified. Let us suppose that (23) holds for j and let us prove it for j + 1. From (8) and from the properties of displacement rank we have

Δ(Â^{(j+1)}) = Δ(Â^{(j)}) − Δ(A_{−1}^{(j)}A_{0}^{(j)^{−1}}A_{1}^{(j)}).  (24)

From (11) it follows that

\begin{matrix}Δ(A_{−1}^{(j)}A_{0}^{(j)^{−1}}A_{1}^{(j)}) = \\ Δ(A_{−1}^{(j)})A_{0}^{(j)^{−1}}A_{1}^{(j)} − A_{−1}^{(j)}A_{0}^{(j)^{−1}}Δ(A_{0}^{(j)})A_{0}^{(j)^{−1}}A_{1}^{(j)} + A_{−1}^{(j)}A_{0}^{(j)^{−1}}Δ(A_{1}^{(j)}). \end{matrix}

By substituting in the above equation relations (16), it may be easily proved, by applying (15), that

\begin{matrix}Δ(A_{−1}^{(j)}A_{0}^{(j)^{−1}}A_{1}^{(j)}) = \\ −v_{−1}^{(j+1)}c_{−1}^{(j)T} A_{0}^{(j)^{−1}}A_{1}^{(j)} − a_{−1}^{(j+1)}u_{0}^{(j+1)T} + A_{−1}^{(j)}A_{0}^{(j)^{−1}}v_{−1}^{(j)}c_{1}^{(j)} + A_{−1}^{(j)}A_{0}^{(j)^{−1}}r^{(j)}u_{0}^{(j)}. \end{matrix}

(25)

Whence, from (24), (25) and from the inductive hypothesis, it follows that

Δ(Â^{(j+1)}) =
−e_{1}c_{−1}^{(0)T} + a_{−1}^{(j+1)}u_{0}^{(j+1)T} − v_{−1}^{(j+1)}e_{m}^{T}(Â^{(j)} − A_{−1}^{(j)}A_{0}^{(j)^{−1}}A_{1}^{(j)})+
(a_{−1}^{(j)} − A_{−1}^{(j)}A_{0}^{(j)^{−1}}r^{(j)})u_{0}^{(j)T} − A_{−1}^{(j)}A_{0}^{(j)^{−1}}v_{−1}^{(j)}c_{1}^{(j)T} .

On the other hand, by following the the same lines of the proof of Theorem 3 and from (17), it may be verified that

(a_{−1}^{(j)} − A_{−1}^{(j)}A_{0}^{(j)^{−1}}r^{(j)})u_{0}^{(j)T} = A_{−1}^{(j)}A_{0}^{(j)^{−1}}v_{−1}^{(j)}c_{1}^{(j)T} .
□

The vectors â^{(j)}, ĉ^{(j)} obtained at two subsequent stages of cyclic reduction are related by means of the equations:

\begin{matrix}â^{(j+1)} = â^{(j)} − A_{−1}^{(j)}A_{0}^{(j)^{−1}}a_{1}^{(j)} \\ ĉ^{(j+1)T} = ĉ^{(j)T} − c_{−1}^{(j)T} A_{0}^{(j)^{−1}}A_{1}^{(j)}, j ≥ 0. \end{matrix}  (26)

Observe that all the above displacement properties still hold when the entries a_{i} of the blocks A_{i}, i = −1, 0, 1, of (1) are square matrices of dimension h > 1. In this case Δ(·) is the block displacement operator, Z is the m × m block down-shift matrix, and the vectors u_{i}^{(j)}, v_{i}^{(j)}, a_{i}^{(j)}, c_{i}^{(j)}, r^{(j)}, â_{i}^{(j)}, ĉ_{i}^{(j)} have m components, each component being a h × h matrix.

In the case of scalar entries a_{i}, the formulae displayed in Theorems 2–5 can be substantially simplified. In fact, from the equation (10), we deduce the following:

Proposition 6 The matrix polynomials A^{(j)}(z) and the vectors r^{(j)}(z), ĉ^{(j)}(z) satisfy the following relations

A^{(j)}(z)^{T} = JA^{(j)}(z)J, ĉ^{(j)}(z) = Jr^{(j)}(z),

where J denotes the permutation (reversion) matrix having unit entries along the anti-diagonal.

Whence, the vectors u^{(j)}(z), v^{(j)}(z), a^{(j)}(z), c^{(j)}(z) are related, as stated by the following:

Proposition 7 For the vectors u^{(j)}(z), v^{(j)}(z), a^{(j)}(z), c^{(j)}(z) we have v^{(j)}(z) = Ju^{(j)}(z), c^{(j)}(z) = Ja^{(j)}(z), thus

Δ(A^{(j)}(z)) = a^{(j)}(z)v^{(j)T} (z)J − v^{(j)}(z)a^{(j)T} (z)J.

In the case of symmetric matrices the formulae can be further simplified, since it can be easily verified that at each step j we have A_{1}^{(j)} = A_{−1}^{(j)T}, A_{0}^{(j)} = A_{0}^{(j)T}, Â^{(j)} = Â^{(j)T}. This simplification allows us to describe the algorithm in terms of the six vectors a_{i}^{(j)}, v_{i}^{(j)}, i = −1, 0, â^{(j)}, r^{(j)}.

The following algorithm enables us to perform a single cyclic reduction step by computing the blocks Â^{(j+1)}, A_{i}^{(j+1)}, i = −1, 0, 1, given the blocks Â^{(j)}, A_{i}^{(j)}, i = −1, 0, 1 by means of few FFT’s of length O(m). For the sake of simplicity we assume that m = 2^{s}, for a positive integer s, and that the entries a_{i} are scalar numbers.

Algorithm.

Input: the vectors r^{(j)}, â^{(j)}, v_{i}^{(j)}, i = −1, 0, a_{i}^{(j)}, i = −1, 0, 1, defining the matrices Â^{(j)}, A_{i}^{(j)}, i = −1, 0, 1, by means of theorems 4,5 and proposition 7;

Output: the vectors r^{(j+1)}, â^{(j+1)}, v_{i}^{(j+1)}, i = −1, 0, a_{i}^{(j+1)}, i = −1, 0, 1, defining the matrices Â^{(j+1)}, A_{i}^{(j+1)}, i = −1, 0, 1, by means of theorems 4,5 and proposition 7;

Computation:

1. Compute the vectors A_{0}^{(j)^{−1}}e_{1}, A_{0}^{(j)^{−1}}a_{−1}^{(j)}, A_{0}^{(j)^{−1}}v_{0}^{(j)} ZJA_{0}^{(j)^{−1}}v_{−1}^{(j)}, that fully

define the 6 vectors representing the matrix A_{0}^{(j)^{−1}}, by solving 4 Toeplitz-like m × m systems.

2. Apply relations (15), (17), (26) in order to compute the vectors r^{(j+1)}, â^{(j+1)}, v_{i}^{(j+1)}, i = −1, 0, a_{i}^{(j+1)}, i = −1, 0, 1. By using the representation formula of A_{0}^{(j)^{−1}} (see theorem 4) and of the matrices A_{1}^{(j)}, A_{−1}^{(j)}, Â^{(j)}, together with the results of section 2.2, this computation can be performed with a finite number of real FFT’s of length 2m.

The above results, integrated with the update of the right hand side b^{(j)} by means of (4), and with a back substitution stage, leads to a direct algorithm for the solution of the finite system T_{n}x = b, where T_{n} is an n × n matrix, with n = mq, q = 2^{p}, p positive integer.

The arithmetic cost of the latter algorithm can be easily estimated. In fact, denoting t(m) the cost of solving an m × m Toeplitz-like system, it follows that O(n log m + (m log m + t(m)) log \frac{n}{m}) arithmetic operations are sufficient to carry out the algorithm. The solution of an m × m Toeplitz-like system can be computed by extending the known algorithms for Toeplitz systems. Therefore, according to the specific properties of the matrices involved (symmetry, positive definiteness, conditioning) algorithms like the fast, super fast, or iterative ones can be used for the computational cost t(m) being O(m^{2}), O(m log^{2} m), O(m log m), respectively.

The cyclic reduction algorithm has also nice properties with respect to parallel computations. In fact, consider the PRAM model where we assume to have several processors working in parallel, at each step a processor performs an arithmetic operation, and the complexity of an algorithm is measured in terms of the number of parallel steps and of processors needed. Then the algorithm can be carried out in O(t_{s}(m) log \frac{n}{m}) parallel steps, using at most max\{t_{p}(m), n\} processors, where t_{s}(m), t_{p}(m) denotes the number of parallel steps and of processors, respectively, needed to solve an m × m Toeplitz-like system. For well conditioned matrices, where the preconditioned conjugate gradient technique can be used, we have t_{s}(m) = log m, t_{p}(m) = m.

4 Convergence properties.

In this section we relate the polynomial a(z) = z^{m} ∑_{k=m_{2}}^{m_{1}} z^{k}a_{k} associated with the matrix T with the polynomial α(z) = det A(z). The results that we obtain are used in the next section to prove some useful convergence properties of the cyclic reduction algorithm.

Let us first recall the following result of [18] that relates the zero of a(z) and α(z).

Theorem 8 For the polynomial α(z) it holds

α(z^{m}) = ∏_{i=0}^{m−1} a(zω^{i}), ω = cos (\frac{2π}{m} ) + i sin (\frac{2π}{m} ) , i^{2} = −1.

A consequence of the above theorem is that, if ξ is zero of the polynomial a(z) then ξ^{m} is zero of α(z). Conversely, if η is zero of α(z) then there exists an m-th root of η which is zero of a(z).

Observe also that both the polynomials a(z) and α(z) have degree m + m_{1}. This is trivial for a(z); in order to prove it for α(z) it is sufficient to consider, in the expansion of det A(z), the term obtained by multiplying the entries on the m_{1}-th diagonal below the main diagonal and the entries in the (m−m_{1})-th diagonal above the main one yielding the z^{m+m_{1}} term.

It is surprising to observe that the cyclic reduction step, is equivalent to the squaring step in the Graeffe algorithm which is customarily used for factoring polynomials [27], [12]. In fact by taking determinants in both sides of (10) we obtain

α^{(j+1)}(z^{2}) = −α^{(j)}(z)α^{(j)}(−z)/ det A_{0}^{(j)}.

An immediate consequence of the above relation is that the polynomials α^{(j)}(z) have degree m + m_{1} and their zeros are explicitly given by

ξ_{i}^{m2^{j}}, i = 1, … m + m_{1},

where

a(ξ_{i}) = 0, i = 1, …, m + m_{1}.

This property is very useful to prove some convergence results of the cyclic reduction algorithm.

First observe that the vector ξ^{T} = (1, ξ, ξ^{2}, …) is such that ξ^{T} T_{∞} = (w^{T} , 0, …), where w^{T} is a suitable vector of m components, for ξ ranging in the set of zeros of a(z). Moreover, assume for simplicity that a(z) has at least m distinct non-null zeros ξ_{1}, …, ξ_{m}, and consider the m × m Vandermonde matrix V_{m} = (ξ_{i}^{j−1})_{i,j=1,m} made up with this set of zeros. It holds

(V_{m}, DV_{m}, D^{2}V_{m}, …)T_{∞} = (W, O, …), D = Diag(ξ_{1}^{m}, …, ξ_{m}^{m}).  (27)

Observe that, for a polynomial a(z) having multiple zeros, the above relation still holds if the matrix V_{m} is replaced by a generalized Vandermonde matrix. The homogeneous part of the above relation can be rewritten as

V_{m}A_{−1} + DV_{m}A_{0} + D^{2}V_{m}A_{1} = O.

By applying the cyclic reduction to the system (27) we deduce that

V_{m}A_{−1}^{(j)} + D^{2^{j}}V_{m}A_{0}^{(j)} + D^{2^{j+1}}V_{m}A_{1}^{(j)} = O, j = 0, 1, …

Multiplying the above relation on the left by V_{m}^{−1} we find

A_{−1}^{(j)} + F^{m2^{j}}A_{0}^{(j)} + F^{m2^{j+1}}A_{1}^{(j)} = O, j = 0, 1, …

where F = V_{m}^{−1}DV_{m} is the Frobenius matrix associated with the polynomial having zeros ξ_{i}, i = 1, …, m.

The following result can be obtained from the above relations.

Theorem 9 Let ξ_{i}, i = 1, …, m + m_{1}, be the zeros of a(z) ordered such that |ξ_{1}| ≤ |ξ_{2}| ≤ … ≤ |ξ_{m+m_{1}}|. If |ξ_{m}| < 1 and if the blocks A_{i}^{(j)} have norm bounded from above by a constant, then lim A_{−1}^{(j)} = O. Moreover, if ξ_{m} is simple then for any operator norm || · || it holds ||A_{−1}^{(j)}|| = O(|ξ_{m}|^{m2^{j}}), otherwise, for any ϵ > 0 there exists an operator norm || · || such that ||A_{−1}^{(j)}|| = O((|ξ_{m}| + ϵ)^{m2^{j}}). Similarly, if |ξ_{m_{1}+1}| > 1 and the blocks have bounded norm then lim A_{1}^{(j)} = O. Moreover it holds ||A_{1}^{(j)}|| = O(|ξ_{m_{1}+1}|^{−m2^{j}}), for any operator norm, or ||A_{1}^{(j)}|| = O((|ξ_{m_{1}+1}|−ϵ)^{−m2^{j}}) for a suitable norm, according to the multiplicity of ξ_{m_{1}+1} as root of a(z).

Proof. For any ϵ > 0 let || · || be an operator norm such that ||F || < ρ(F) + ϵ (compare [19]), where ρ(F) denotes the spectral radius of F.. Then we have

||A_{−1}^{(j)}|| ≤ ||F ||^{m2^{j}}||A_{0}^{(j)}|| + ||F ||^{m2^{j+1}}||A_{1}^{(j)}||,

whence, since ||A_{0}^{(j)}||, ||A_{1}^{(j)}|| are bounded from above by a constant, it holds

||A_{−1}^{(j)}|| = O((ρ(F) + ϵ)^{m2^{j}}) = O((|ξ_{m}| + ϵ)^{m2^{j}}).

If ξ_{m} is a simple root of a(z) then there exists a norm || · || such that ||F || = |ξ_{m}|. Therefore, by applying the same argument as before the above relation turns into ||A_{−1}^{(j)}|| = O(|ξ_{m}|^{m2^{j}}). In order to prove that ||A_{1}^{(j)}|| = O((|ξ_{m_{1}+1}| − ϵ)^{−m2^{j}}) we apply a similar argument to the relation A_{−1}^{(j)} + F̂^{m2^{j}}A_{0}^{(j)} + F̂^{m2^{j+1}}A_{1}^{(j)} = O, where F̂ is the Frobenius matrix associated with the polynomial having zeros ξ_{m_{1}+1}, …, ξ_{m_{1}+m}. □

A nice consequence of theorem 9 is that, if T^{(0)} = T_{∞} or T^{(0)} = T_{±∞}, the sequence T^{(j)} converges to a block (bi-)diagonal matrix. Thus, for j = O(log log ϵ^{−1}), we may approximate an arbitrary finite number of components of the solution of the system T^{(j)}x^{(j)} = b^{(j)}, within the error ϵ. If our algorithm is applied to solve the finite system T_{n}x = b, where n = m2^{p} and p >> log log ϵ^{−1}, it is sufficient to perform j = O(log log ϵ^{−1}) steps of cyclic reduction, and then to compute the solution of T^{(j)}x^{(j)} = b^{(j)}, within the error ϵ, by solving a block (bi-)diagonal system. In this case the computational cost is O(n log m + t(m) log log ϵ^{−1}).

Conditions under which the blocks A_{i}^{(j)} have bounded norm are given in [22], [33]. These conditions can be substantially relaxed in the real symmetric case where we only require that the polynomial a(z) has no zeros of modulus 1. In this case the convergence speed depends on the modulus of the zero ξ closest to the unit circle. These properties are investigated with more details in the next section.

A very important situation where the blocks A_{i}^{(j)} have bounded norm is the case T_{n} = I − P_{n}, with P_{n} stochastic matrix. This case occurs in the analysis of Markov chains arising in queueing problems [9], [10].

4.1 The real symmetric positive definite case.

Consider the case where the matrix T_{n} is real symmetric, so that for the blocks A_{i}^{(j)} it holds A_{−1}^{(j)} = A_{1}^{(j)T} and A_{0}^{(j)} = A_{0}^{(j)T}, j = 0, 1, . . .. Assume for simplicity that the polynomial a(z) has no zero of modulus 1. Due to the symmetricity properties, the polynomial a(z) is symmetric, that is, a(z) = z^{2m}a(z^{−1}), and its zeros occur in pairs ξ_{i} = ξ_{2m−i+1}^{−1}, i = 1, …, m. Let us arrange the zeros ξ_{i} in such a way that |ξ_{1}| ≤ |ξ_{2}| ≤ … ≤ |ξ_{m}| < 1 < |ξ_{m+1}| ≤ … ≤ |ξ_{2m}|.

Due to the Szëgo theory on Toeplitz matrices [21], the matrices T_{n}, are either positive definite or negative definite for any n according to the sign that a(z) takes over the unit circle. Throughout we assume, for simplicity T_{n} positive definite. Moreover, the eigenvalues of the matrix T_{n} belong to the interval [min_{|z|=1} |a(z)|, max_{|z|=1} |a(z)|].

It is useful to observe that the block tridiagonal matrices T^{(j)}, generated at each step of the cyclic reduction, being Schur complements in T_{n}, have their eigenvalues belonging to the interval [min λ(T_{n}), max λ(T_{n})], where λ(A) denotes the set of

eigenvalues of the matrix A. Similarly the diagonal blocks A_{0}^{(j)}, Â^{(j)}, [ \begin{matrix}A_{0}^{(j)} & A_{1}^{(j)T} \\ A_{1}^{(j)} & A_{0}^{(j)} \end{matrix} ],

share the same properties due to the Cauchy interlace theorem [19]. Therefore it follows that the all the blocks A_{0}^{(j)}, A_{1}^{(j)}, Â^{(j)}, have norm bounded from above by a constant.

These observations allow us to conclude with the following theorem

Theorem 10 Let T_{n} be a real symmetric matrix (Hermitian) matrix such that the associated polynomial a(z) has no zeros in the unit circle. Then for the matrices A_{0}^{(j)}, A_{1}^{(j)}, Â^{(j)}, generated by the cyclic reduction we have

lim A_{1}^{(j)} = O, ||A_{1}^{(j)}|| ≤ η|ξ_{m}|^{m2^{j}}

The above theorem can be used in order to devise an algorithm for the approximation of n components of the solution of the system T x = b, in O((m log m + t(m)) log log ϵ^{−1}+n log m) arithmetic operations provided that n is less than the size of T. More precisely, in the light of theorem 9, we may replace the term log log ϵ^{−1} with the expression log_{2} log θϵ^{−1} − log_{2} log |ξ_{m−1}|^{−m} for a suitable constant θ depending on the specific norm adopted in the proof of theorem 9.

Another consequence of the boundness of the eigenvalues of T^{(j)} is that cond(T^{(j)}) ≤ cond(T) where cond is the condition number expressed in terms of the Euclidean norm. This means that the intermediate systems that are encountered in the cyclic reduction steps are not worstly conditioned than the initial system.

4.2 The unsymmetric/indefinite case: normal equations

If the polynomial a(z) is not symmetric or if a(ξ) = 0 for some ξ of unit modulus, the convergence results of the previous section do not generally hold and the cyclic reduction technique cannot be used for solving the linear system T_{n}x = b in an iterative fashion. Moreover, we cannot bound with a constant the condition number of the m × m blocks to be inverted at each step of cyclic reduction unless a(z) is positive for |z| = 1. Therefore, the stability of the cyclic reduction technique is not insured in this case.

However, with a simple trick we may reduce the general case to the positive definite one.

Let us first assume that a(z) has zeros ξ_{1}, … ξ_{m+m_{1}} and that |ξ_{k}| = 1 for a given k. We may replace the system T_{n}x = b with the systems

\begin{matrix}DT_{n}D^{−1}y = Db \\ y = Dx, \end{matrix}

where D = Diag(1, γ, γ^{2}, …, γ^{n−1}), for a suitable γ. In this way the new matrix DT_{n}D^{−1} is still Toeplitz and is associated with the polynomial a_{γ}(z) = a(γz). The parameter γ can be chosen in such a way that a_{γ}(z) has no zero of modulus 1. This can be done, deterministically, if informations about the moduli of the zeros of a(z) are known, or by means of randomization with a random choice of γ in the range [1, 1 + σ] for a given positive σ, which guarantees with probability 1 that a_{γ}(z) has no zero on the unit circle.

Indeed, if T_{n} is real symmetric, the new system is not anymore symmetric. In order to keep the symmetricity property, we switch to the normal equations

\begin{matrix}K_{n}y = c \\ K_{n} = (DT_{n}D^{−1})^{T} (DT_{n}D^{−1}) \\ c = (DT_{n}D^{−1})^{T} Db. \end{matrix}  (28)

It is easy to check that K_{n} has 4m + 1 diagonals and is still Toeplitz except for the entries in the leading and in the trailing m × m sub-matrices. Moreover, for the Toeplitz part, the matrix is associated with the polynomial q_{γ}(z) = ϕ_{γ}(z)ϕ_{γ}(z^{−1})z^{(m_{1}+m)} having zeros ξ_{1}, …, ξ_{m+m_{1}}, ξ_{1}^{−1}, …, ξ_{m_{1}+m}^{−1}. More precisely K_{n} = V_{n} − Diag(R, O, S), where V_{n} is the banded Toeplitz matrix associated with q_{γ}(z), R = A_{−1}^{T}A_{0}^{−1}, S = A_{1}^{T} A_{1}, and O denotes a null matrix of size n − 2m.

Partitioning K_{n} into 2m × 2m blocks yields

K_{n} = [ \begin{matrix}Ĉ & B_{−1} \\ B_{1} & B_{0} & B_{−1} \\ & ... & ... & ... \\ & & B_{1} & B_{0} & B_{−1} \\ & & & B_{1} & C̃ \end{matrix} ]  (29)

Let us assume n = 2mq, where q = 2^{p} − 1 and p is a positive integer. In this way the cyclic reduction algorithm applied to the system (28) generates a sequence of systems whose matrices, of dimension (2^{p−j} − 1)2m, have the same structure of (29) and are defined by the 2m × 2m blocks Ĉ^{(j)}, C̃^{(j)}, B_{i}^{(j)}, i = −1, 0, 1, such that

\{ \begin{matrix}B_{−1}^{(j+1)} = −B_{−1}^{(j)}B_{0}^{(j)^{−1}}B_{−1}^{(j)} \\ B_{0}^{(j+1)} = B_{0}^{(j)} − B_{1}^{(j)}B_{0}^{(j)^{−1}}B_{−1}^{(j)} − B_{−1}^{(j)}B_{0}^{(j)^{−1}}B_{1}^{(j)} \\ B_{1}^{(j+1)} = −B_{1}^{(j)}B_{0}^{(j)^{−1}}B_{1}^{(j)} \\ B̂^{(j+1)} = B̂^{(j)} − B_{−1}^{(j)}B_{0}^{(j)^{−1}}B_{1}^{(j)} \\ Ĉ^{(j+1)} = B̂^{(j+1)} + Ĉ − B_{0} \\ C̃^{(j+1)} = JB̂^{(j+1)}J + C̃ − B_{0} \end{matrix}

In the light of the results of the previous section the matrices Ĉ^{(j)}, C̃^{(j)} and B_{i}^{(j)}, i = −1, 0, 1, have displacement rank independent of j and an algorithm similar to the one of section 4.1 applies.

Observe that the condition number of B_{0}^{(j)} is bounded by the condition number of V_{n} and thus depends on max_{|z|=1} q_{γ}(z)/ min_{|z|=1} q_{γ}(z). A similar bound cannot be proved for the blocks Ĉ^{(j)},C̃^{(j)}.

5 Numerical experiments and applications.

We have implemented several versions of our algorithms in Fortran 90. The most general one deals with banded (unsymmetric) block Toeplitz matrices. Two more specific versions deal with the cases of (unsymmetric) scalar Toeplitz matrices and with symmetric scalar Toeplitz matrices. For the inversion of the diagonal blocks we used the conjugate gradient method with no preconditioning. We performed several numerical experiments, here we report the ones concerning the scalar symmetric case.

The test matrices that we have considered are defined by the function ϕ(z) = ∑_{i=−m}^{m} z^{i} + (a_{0} − 1), for m = 16, 32, 64, n = 2^{h}, h = 14, 15, …, 22, where the value of a_{0} has been chosen differently in such a way that the polynomial z^{m}ϕ(z) has some zeros more or less close to the unit circle. The right-hand vector of the system has been chosen so that the system has the solution x = (1, …, 1)^{T}.

More specifically, for m = 16, we tested the algorithm with a_{0} = 8.2. In this case, the polynomial ϕ(z)z^{m} has no zeros in the unit circle, the zeros closest to the unit circle are two pairs of modulus 1.003052 and 0.99695731. Their ratio, 1.006, is far enough from 1 in order to allow the double exponential convergence of the cyclic reduction algorithm.

For m = 32 we have chosen a_{0} = 15.1315. With this value, the polynomial ϕ(z)z^{m} has no zeros of modulus 1. The zeros closest to the unit circle have modulus 1.00009 and 0.9999099, and their ratio has modulus 1.00018.

For m = 64 we set a_{0} = 25. With this value the cyclic reduction algorithm does not converge since the polynomial ϕ(z)z^{m} has two pairs of complex zeros of modulus 1.

CR  PCG
n  Time (s.) Residual Time (s.) Residual PCG/CR
8192  0.6  3. e -13  1.0  2. e -11  1.7
16384  1.0  3. e -13  2.6  2. e -11  2.6
32768  1.8  7. e -13  6.9  2. e -11  3.8
65536  3.5  7. e -13  16.0  2. e -11  4.6
131072  6.7  7. e -13  36.9  2. e -11  5.5
262144  13.3  7. e -13  132.0  2. e -11  9.9
524288  26.5  7. e -13  376.3  2. e -11  14.2
1048576  52.8  7. e -13  *  *  *
2097152  98.3  7. e -13  *  *  *
Table 1: m = 16, a_{0} = 8.2

CR  PCG
n  Time (s.) Residual Time (s.) Residual PCG/CR
8192  0.7  2. e -11  1.3  1. e -12  1.9
16384  1.1  1. e -10  3.3  1. e -12  3.0
32768  1.8  2. e -10  7.9  1. e -12  4.3
65536  3.4  2. e -10  17.2  4. e -11  5.1
131072  6.4  2. e -10  38.9  4. e -11  6.1
262144  12.3  2. e -10  139.3  4. e -11  11.3
524288  24.3  2. e -10  395.9  4. e -11  16.3
1048576  48.2  2. e -10  *  *  *
2097152  95.6  2. e -10  *  *  *
Table 2: m = 32, a_{0} = 15.1315

The algorithm has been compared with the preconditioned conjugate gradient method obtained by using the circulant preconditioner of G. Strang [32], on an alpha workstation.

In the tables 1, 2, 3 we report the time in seconds needed by the two algorithms (a “*” denotes failure of the algorithm due to lack of memory), the maximum error in the solution and the ratios between the times needed.

The numerical experiments confirm that our algorithm is faster than the preconditioned conjugate gradient method, more efficient in terms of memory use and numerically stable.

References

[1] G. S. Ammar and B. Gragg. Numerical experience with superfast real Toeplitz solver. Linear Algebra Appl., 121:185–206, 1989.

[2] O. Axelsson and G. Lindskög. The rate of convergence of the preconditioned conjugate gradient method. Numer. Math., 52:499–523, 1986.

[3] D. Bini. On a class of matrices related to Toeplitz matrices. Technical Report 83-5, SUNY University, Albany, New York, 1983.

[4] D. Bini. Parallel solution of certain Toeplitz linear systems. SIAM J. Comput., 13:268–276, 1984.

CR  PCG
n  Time (s.) Residual Time (s.) Residual PCG/CR
8192  0.8  2. e -11  1.0  3. e -12  1.3
16384  1.2  4. e -11  2.5  1. e -11  2.1
32768  2.1  1. e -11  6.6  2. e -11  3.1
65536  3.6  6. e -12  15.4  6. e -12  4.3
131072  6.8  2. e -10  37.0  1. e -11  5.4
262144  13.0  4. e -11  123.0  4. e -11  9.5
524288  25.4  4. e -11  328.0  6. e -12  12.9
1048576  49.8  1. e -9  *  *  *
2097152  99.2  8. e -11  *  *  *
Table 3: m = 64, a_{0} = 25

[5] D. Bini. Matrix structures in parallel matrix computations. Calcolo, 25:37–51, 1988.

[6] D. Bini and M. Capovani. Spectral and computational properties of band symmetric Toeplitz matrices. Linear Algebra Appl., 52:99–126, 1983.

[7] D. Bini and F. Di Benedetto. A new preconditioner for the parallel solution of positive definite Toeplitz systems. In Proceedings 2^{nd} Annual SPAA, Crete, Greece, pages 220–223. ACM Press, July 1990.

[8] D. Bini and B. Meini. On cyclic reduction applied to a class of Toeplitz-like matrices arising in queueing problems. In W. J. Stewart, editor, Computations with Markov Chains, pages 21–38. Kluwer Academic Publisher, 1995.

[9] D. Bini and B. Meini. On the solution of a nonlinear matrix equation arising in queueing problems. SIAM J. Matrix Anal. Appl., 17:906–926, 1996.

[10] D. Bini and B. Meini. Improved cyclic reduction for solving queueing problems. Numerical Algorithms, 15:57–74, 1997.

[11] D. Bini and B. Meini. Inverting block Toeplitz matrices in block Hessenberg form by means of displacement operators: application to queueing problems. Linear Algebra Appl., 272:1–16, 1998.

[12] D. Bini and V. Pan. Matrix and Polynomial Computations, Vol. 1: Fundamental Algorithms. Birkhäuser, Boston, 1994.

[13] S. Bondeli and W. Gander. Cyclic reduction for special tridiagonal systems. SIAM J. Matrix Anal. Appl., 15:321–330, 1994.

[14] R. H. Chan. Circulant preconditioners for Hermitian Toeplitz systems. SIAM J. Matrix Anal. Appl., 10:542–550, 1989.

[15] R. H. Chan and M. K. Ng. Conjugate gradient methods for Toeplitz systems. SIAM Review, 38:427–482, 1996.

[16] R. H. Chan and M. K. Ng. Scientific applications of iterative Toeplitz solvers. Calcolo, 33:249–267, 1996.

[17] R. H. Chan and P. Tang. Fast band-Toeplitz preconditioners for Hermitian Toeplitz systems. SIAM J. Sci. Comput., 15:164–171, 1994.

[18] H. R. Gail, S. L. Hantler, and B. A. Taylor. Non-skip-free M/G/1 and G/M/1 type Markov chains. Adv. Appl. Prob., 29:733–758, 1997.

[19] G.H. Golub and C.F. van Loan. Matrix Computations. The Johns Hopkins University Press, Baltimore, 1989.

[20] J. Grcar and A. Sameh. On certain parallel Toeplitz linear system solvers. SIAM J. Sci. Stat. Comput., 2:238–256, 1981.

[21] U. Grenander and G. Szegö. Toeplitz Forms and Their Applications. Chelsea Publishing, New York, 2nd edition, 1984.

[22] Don Heller. Some aspects of the cyclic reduction algorithm for block tridiagonal linear systems. SIAM J. Num. Anal., 13:484–496, 1976.

[23] T. Kailath and A.H. Sayed. Displacement structure: theory and applications. SIAM Review, 37:297–386, 1995.

[24] T. Kailath, A. Viera, and M. Morf. Inverses of Toeplitz operators, innovations and orthogonal polynomials. SIAM Review, 20:106–119, 1978.

[25] M. F. Neuts. Matrix-Geometric Solutions in Stochastic Models. Johns Hopkins University Press, Baltimore, 1981.

[26] M. F. Neuts. Structured Stochastic Matrices of M/G/1 Type and Their Applications. Marcel Dekk., New York, 1989.

[27] M. A. Ostrowski. Recherches sur la methode de Graeffe et les zeros des polynomes et des series de Laurent. Acta Math., 72:99–257, 1940.

[28] S. Serra. Preconditioning strategies for asymptotically ill-conditioned block Toeplitz systems. BIT, 34:579–594, 1994.

[29] S. Serra. Optimal, quasi-optimal and superlinear preconditioners for asymptotically ill-conditioned positive definite Toeplitz systems. Math. Comp., 66:651– 665, 1997.

[30] S. Serra Capizzano. Toeplitz preconditioner constructed from linear approximation processes. in printing, 1998.

[31] J. Stoer and R. Bulirsch. Introduction to Numerical Analysis. Springer Verlag, New York, 1980.

[32] G. Strang. A proposal for Toeplitz matrix computations. Studies in Applied Mathematics, 74:171–176, 1986.

[33] P. Yalamov and V. Pavlov. Stability of the block cyclic reduction. Linear Algebra Appl., 249:341–358, 1996.