SOLVING BLOCK BANDED BLOCK TOEPLITZ SYSTEMS WITH STRUCTURED BLOCKS: ALGORITHMS AND APPLICATIONS

DARIO ANDREA BINI AND BEATRICE MEINI^{∗}

Abstract. After motivating the interest and the role of Toeplitz matrices by means of several examples of applications, we focus our interest on the solution of block banded block Toeplitz systems endowed with a two-level structure. In particular we consider block tridiagonal block Toeplitz systems where the blocks are Toeplitz themselves or have a Kronecker product structure, or are Frobenius-like matrices. We introduce some computational tools like the orthogonal displacement representation and the concept of approximate displacement rank and combine them with other well known tools for the design of effective algorithms for the solution of structured problems. Applications to polynomial factorization, solution of queueing models, and image restoration are shown. Finally the concept of approximate displacement rank is used for the design of effective algorithms for Toeplitz-like matrix inversion based on Newton’s iteration.

1. Matrix structures and their role. The mathematical modeling of problems of the real world often leads to solving linear systems of equations endowed with specific structures and with very large dimensions. Generally the matrices involved in these models are defined by few parameters that characterize all their entries through a given pattern.

Certain structures are very frequent and reflect specific features that are common to problems arising in diverse fields of applied mathematics and engineering. Toeplitz matrices are an example of this fact. A matrix is Toeplitz if it has constant entries along each diagonal. This kind of matrices, with their block versions, arises whenever properties of shift invariance are satisfied by some function in the model. They are encountered, in particular, in fields like image processing, queueing theory, computer algebra and in the numerical solution of differential equations where the shift invariance takes different forms.

Indeed, the stronger is the structure, the more likely efficient solution algorithms can be designed. A challenging problem is to analyze and exploit, as most as possible, the specific structures that characterize a problem in order to devise tools for the design of efficient solution algorithms.

Very often such structures are accompanied by further properties, like banded patterns or Kronecker product patterns that make the problem still more interesting. For block matrices we may have structures at two (or more) different levels: an inner structure, typical of each block, and an outer structure proper of the block matrix. Exploiting both structures is a nontrivial task as it could seem at first.

In fact, it may happen that algorithms designed for the exploitation of the outer structure destroy the inner one and vice-versa.

Concerning Toeplitz computation there is a wide and consolidated literature that provides a lot of tools for the design of efficient algorithms.

In this paper we first describe some examples of applications of (block) Toeplitz matrices and relate the specific structures encountered in this analysis with the peculiar features of the problem. Then we present some recent algorithms, based on the Toeplitz matrix technology, for the solution of certain structured linear systems. In particular, we consider matrices with a two-level structure, like block banded block Toeplitz matrices with banded Toeplitz blocks, or with blocks having a Kronecker

∗ Dipartimento di Matematica, Università di Pisa, via Buonarroti 2, I-56127 Pisa, Italy. E-mail:

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

product pattern. These problems arise in queueing theory, in image processing and in polynomial computations.

We provide some new techniques for dealing with such problems and discuss some related issues. In particular we introduce the concepts of orthogonal displacement representation, useful to improve the numerical stability of the computation, and the concept of approximate displacement rank. The latter tool is particularly suitable for matrices having a small Euclidean distance from the set of matrices with low displacement rank.

We show applications of these concepts to improving and extending known techniques, like cyclic reduction, to multilevel Toeplitz matrices. We also apply approximate displacement rank for the acceleration of Newton’s iteration for inverting Toeplitz-like matrices.

The paper is organized in the following way. This section continues with some examples where the Toeplitz structure is fundamental in the modeling of the mathematical problem. Section 2 is devoted to the description of some tools for structure analysis, in particular we introduce orthogonal displacement representations, approximate displacement rank and recall the main features of cyclic reduction.

In section 3 we present the solution of diverse problems; specifically, the solution of banded Toeplitz systems with application to polynomial factorization; the solution of two queueing models; the solution of block banded block Toeplitz systems encountered in image processing; the inversion of Toeplitz-like matrices by means of Newton’s iteration and approximate displacement rank.

1.1. Structures in queueing problems. In the numerical solution of Markov chains, we need to solve (infinite) linear systems Qx = b, where Q = I − P and P is a stochastic, or substochastic matrix. The matrix P, whose (i, j) entry p_{i,j} is the probability of transition from the state i to the state j of the system in one unit of time, can be strongly structured, according to the specific model.

In many problems of the real world (eg. network communication, queueing problems), the transition probabilities p_{i,j} are shift invariant, i.e., p_{i+h,j+h} = p_{i,j}. This is a very natural situation in queueing models where the probability that the number j of customers waiting in the system is changed into i, in a time unit, depends only on the difference i − j, that coincides with the difference between the number of new arrivals and the number of customers served in the unit time. Evidently this shift invariance property leads to Toeplitz structured matrices.

In many models it is not allowed that in a time unit the number of waiting customers is increased by an additive number larger than k or smaller than −h, for positive h, k. This property, i.e., p_{i,j} = 0 for j − i > k, or j − i < −h, leads to a lower, or upper, generalized Hessenberg structure.

For systems with more than one queue, the scalar Toeplitz, lower/upper (generalized) Hessenberg structure of the transition probability matrix is replaced with a block Toeplitz, block lower/upper (generalized) Hessenberg structure.

In particular, in the case of QBD (Quasi-Birth-Death) processes [26], Q = I − P is a block tridiagonal block Toeplitz matrix except for the block in position (1,1), i.e.,

Q = [ \begin{matrix}Â & C & & ◯ \\ B & A & C \\ & B & A & ... \\ ◯ & & ... & ... \end{matrix} ] , (1.1)

where A, Â, B, C are m × m matrices.

The block entries A, Â, B, C can have also an inner structure; for example, for PH/PH/1 queues [26, 3], that are a particular case of QBD processes, the blocks, having size m = m_{1}m_{2}, are given by

A = T ⊕ S, C = T^{0}α ⊗ I_{m_{2}}, B = I_{m_{1}} ⊗ S^{0}β, (1.2)

where T ∈ R^{m_{1}×m_{1}}, S ∈ R^{m_{2}×m_{2}}, α ∈ R^{1×m_{1}}, β ∈ R^{1×m_{2}}, T^{0} and S^{0} are such that T e + T^{0} = 0 and Se + S^{0} = 0, e = (1, 1, …, 1)^{T}, the symbols ⊗ and ⊕ denote the Kronecker product and the Kronecker sum of matrices, respectively, and I_{r} denotes the r × r identity matrix.

Another interesting class of queueing problems is the one modeled by Non-skip-free M/G/1-type Markov chains [16, 8]; in this case the matrix Q is a block Toeplitz matrix in generalized block Hessenberg form, i.e.,

Q = [ \begin{matrix}A_{k} & A_{k−1} & … & A_{0} & & ◯ \\ A_{k+1} & A_{k} & A_{k−1} & … & A_{0} \\ A_{k+2} & A_{k+1} & A_{k} & ... & ... & ... \\ ⋮ & ... & ... & ... & ... & ... \end{matrix} ] , (1.3)

where A_{i}, i = 0, 1, …, are m × m matrices. In the applications, the above matrix is approximated with a banded matrix, i.e., it is assumed that A_{i} = 0, for i > N, for a suitable integer N. The matrix Q can be reblocked into k × k block matrices thus obtaining a block Hessenberg block Toeplitz matrix with structured block entries. Hence, solving Qx = b can be reduced to the problem of the solution of a system with a block Toeplitz block Hessenberg matrix with structured blocks.

Infinite block Toeplitz block Hessenberg systems can be reduced to block tridiagonal block Toeplitz systems where the matrix, of the kind (1.1), has infinite blocks, by means of a procedure devised by V. Ramaswami in [28] that consists in adding new auxiliary variables and new equations. The infinite blocks A, B, C that are obtained in this way have a strong structure that is related to block Frobenius matrices. More precisely, if the original block Hessenberg block Toeplitz matrix is defined by (1.3) with k = 1, then we have

\begin{matrix}A = [ \begin{matrix}A_{1} & A_{2} & A_{3} & … \\ & I \\ & & I \\ . & . & . \end{matrix} ] , B = [ \begin{matrix}A_{0} & 0 & … \\ 0 & 0 & … \\ ⋮ & ⋮ \end{matrix} ] , \\ C = [ \begin{matrix}0 & 0 & … \\ −I & 0 & ... \\ & −I & ... \\ . & . & . \end{matrix} ] . \end{matrix}

Once again we encounter a matrix with a two-level structure where the external block Toeplitz structure is combined with the internal block Frobenius structure.

1.2. Structures in polynomial computations. Polynomials, power series and (infinite) Toeplitz matrices are closely related. If a(x) = ∑_{i=0}^{+∞} a_{i}x^{i}, b(x) = ∑_{i=0}^{+∞} b_{i}x^{i}, c(x) = ∑_{i=0}^{+∞} c_{i}x^{i} are formal power series such that c(x) = a(x)b(x), then the relations that express the coefficient of c(x) with those of a(x) and b(x) are synthesized by

[ c_{0} ; c_{1} ; ⋮ ] = [ \begin{matrix}a_{0} \\ a_{1} & a_{0} \\ ⋮ & ... & ... \end{matrix} ] [ b_{0} ; b_{1} ; ⋮ ] .

The same relation holds for polynomials assuming that the coefficients of x^{i} in the power series a(x), b(x) and c(x) are zero when i is larger than the degree. The Toeplitz structure of the matrix in the above relation comes from the fact that multiplying a power series by the variable x corresponds to multiplying its vector coefficients by the Toeplitz matrix Z = (z_{i,j}), where z_{i+1,i} = 1, z_{i,j} = 0 elsewhere.

For matrix power series, i.e., power series where the coefficients are m×m matrices the above relation still holds where the matrix involved turns into a block Toeplitz matrix.

This correlation between Toeplitz matrices and power series allows one to express almost any polynomial or power series computation (like gcd, lcm, Euclidean scheme, power series arithmetic, etc) in terms of Toeplitz computations and vice-versa (see [13], [7], [2]).

Recently, this correlation has been used for devising an efficient algorithm for polynomial and power series factorization [4].

Let p(z) = ∑_{i=0}^{2m} p_{i}z^{i} be a polynomial of degree (at most) 2m, having zeros z_{i}, i = 1, …, 2m, such that |z_{1}| ≤ … ≤ |z_{m}| < |z_{m+1}| ≤ … ≤ |z_{2m}|, and denote f^{T} = (f_{0}, …, f_{m−1}) the coefficient vector of the polynomial f(z) = ∏_{i=1}^{m}(z − z_{i}) = ∑_{i=0}^{m} f_{i}z^{i}, f_{m} = 1.

Consider the n × n banded Toeplitz matrix Q_{n} = (q_{i,j}), q_{i,j} = p_{i−j+m} if 0 ≤ i − j + m ≤ 2m, q_{i,j} = 0, otherwise. A consequence of Koenig’s theorem [19] implies that for n large enough Q_{n} is nonsingular, moreover, denoting w^{(n)} the vector made up by the first m components of the first column of Q_{n}^{−1}, it holds f = −Uw^{(n)} +O(σ^{mn}), where σ = |z_{m}/z_{m+1}| and U is the m × m upper triangular Toeplitz matrix defined by its first row p_{2m}, …, p_{m+1} (compare with [4]).

Thus the problem of approximating a factor of a polynomial is reduced to compute the first m entries of the first column of the inverse of Q_{n} for a large enough value of n. This computation can be performed by an algorithm that provides the above values for n = m2^{i}, i = 1, 2, 3, …, and where the cost of switching from n = m2^{i} to n = m2^{i+1} is O(m log^{2} m) arithmetic operations.

1.3. Structures in image processing. In many problems of image restoration one of the main issues is to recover an image from a blurred one, where the blur operator is linear, space (shift) invariant and has a local action. More precisely, if we represent images as m × n arrays X of reals, each of them representing a pixel (picture element) of the image, the discrete blur operator can be described as B = L(X), where the blurred image B = (b_{i,j})_{i=1,m,j=1,n} is related to the original image X = (x_{i,j})_{i=1,m,j=1,n} by

b_{i,j} = ∑_{h=−M}^{M} ∑_{k=−M}^{M} α_{h,k}x_{i+h,j+k}, (1.4)

(here we assume that x_{i,j} = 0 if i or j is out of range) and the Point Spread Function (PSF) matrix α = (α_{i,j})_{i,j=−M,...,M} is such that ∑_{i,j} α_{i,j} = 1, α_{i,j} ≥ 0.

The entry α_{i−h,j−k} represents the influence in the pixel (i, j) of the blurred image of a unit source placed in the pixel (h, k) of the original image.

Rearranging the entries of X and B column-wise as vectors x and b in (1.4) we obtain the linear system Qx = b, where the matrix Q is an n × n block Toeplitz block banded matrix with m×m banded Toeplitz blocks. More precisely, Q is a block (2M + 1)-diagonal matrix with (2M + 1)-diagonal blocks. Moreover the blocks on the i-th block diagonal, i = −M, …, M, coincide with the (2M + 1)-diagonal Toeplitz matrix having on the j-th diagonal the entry α_{i,j}, for j = −M, …, M. The block Toeplitz structure is derived by the shift invariance properties of the PSF, i.e., the mask α is independent of the specific pixel to which is applied. The block banded structure is due to the locality of the PSF.

1.4. Structures in differential equations. Differential equations with boundary conditions can be efficiently treated by means of numerical techniques. For instance, the Sturm-Liouville problem

\{ \begin{matrix}−(a(x)u′(x))′ = f(x) \\ u(0) = α_{0}, u(1) = α_{1} \end{matrix}

can be numerically solved by replacing the unknown function u(x) with the unknown vector u = (u_{i}), where u_{i} = u(x_{i}) and x_{i} = i/(n + 1), i = 1, …, n. The differential operator, discretized by means of the finite differences technique, is reduced to a tridiagonal matrix T = (t_{i,j}), t_{i,i} = (a_{i} + a_{i−1})/(n + 1), t_{i,i+1} = t_{i+1,i} = −a_{i}/(n + 1) and the vector u solves the system T u = f + O(h^{2}), h = 1/(n + 1), where f = (f_{i}), f_{i} = f(x_{i}), i = 2, …, n − 1, f_{1} = f(x_{1}) + (n + 1)^{2}α_{0}, f_{n} = f(x_{n}) + (n + 1)^{2}α_{1}.

It is immediate to verify that if a(x) is constant, that is if the differential operator reduces to the second derivative, then the matrix T is Toeplitz; for a generic a(x) the matrix T is not Toeplitz but is “locally Toeplitz” in the sense of [29]. For derivatives of higher order with constant coefficients the matrix T is still banded Toeplitz and its band-width is proportional to the order of the derivative. For a linear differential operator of functions in m variables with constant coefficients, like the Laplace operator, the matrix that provides the numerical discretization over [0, 1]^{m} is a multilevel block Toeplitz matrix with a multilevel band.

The band Toeplitz structure of the matrices discretizing linear differential operators with constant coefficients is essentially due to the shift invariance property of derivatives expressed by the following relation

D^{k}(f(x + h)) = g(x + h), g(x) = D^{k}(f(x)),

valid for the differential operator D(f(x)) = f′(x).

2. Some tools for structure analysis. In this section we recall some basic tools for the design of efficient algorithms for solving certain structured systems. Throughout this section we assume that all the matrices have real entries unless explicitly specified.

2.1. Displacement rank. The concept of displacement operators and displacement rank, introduced in [22, 21] and elaborated in [17, 13], is a powerful tool for dealing with Toeplitz matrices. Here we recall the main ideas and results.

Let us consider the operator

(2.1) Δ(A) = AZ − ZA, Z = (z_{i,j}), z_{i+1,i} = 1, z_{i,j} = 0, otherwise, applied to n × n matrices A, and define the displacement rank of A, associated with Δ, as Drk(A) = rank(Δ(A)). A Toeplitz matrix has displacement rank at most 2, moreover, Drk(A) is a sort of measure of the distance of A from the algebra of lower triangular Toeplitz matrices, for which the displacement rank is zero.

Define L(a) the n × n lower triangular Toeplitz matrix having a as first column. It holds [13]

Δ(A) = ∑_{i=1}^{r} u_{i}v_{i}^{T} ⇔ A = L(Ae_{1}) + ∑_{i=1}^{r} L(u_{i})L^{T} (Zv_{i}), (2.2)

where u_{i}, v_{i}, i = 1, …, r, are n-dimensional vectors, e_{1} denotes the first column of the identity matrix. The relation in the right-hand side of (2.2) is called the displacement representation of A and the vectors a = Ae_{1}, u_{i}, v_{i}, i = 1, …, r are called a set of generators for the displacement representation of A. Equivalent representations can be given in terms of different operators, see for instance [20, 13].

A simple consequence of the relation

Δ(A^{−1}) = −A^{−1}Δ(A)A^{−1}, (2.3)

that can be directly derived from (2.1), is that Drk(A^{−1}) = Drk(A) and the following displacement representation of A^{−1} holds

A^{−1} = L(A^{−1}e_{1}) − ∑_{i=1}^{r} L(A^{−1}u_{i})L^{T} (ZA^{−T} v_{i}). (2.4)

The representations (2.2), (2.4) allow us to make elementary operations among Toeplitz-like matrices at a low cost by means of FFT. In particular, computing the product of a matrix with displacement rank r and a vector costs O(rn log n) arithmetic operations (ops), inverting a strongly nonsingular matrix with displacement rank r costs O(r^{2}n log^{2} n) ops by using a superfast algorithm like the Bitmead Anderson algorithm [14].

It is a simple matter to show that Δ(AB) = AΔ(B) + Δ(A)B, therefore it holds Drk(AB) ≤ Drk(A) + Drk(B). Similarly we have Drk(A + B) ≤ Drk(A) + Drk(B).

2.2. Orthogonal displacement representation. It is worth to point out that the set of generators for the displacement representation of A is not unique. Indeed, observe that, if Δ(A) = UV^{T}, for U, V ∈ R^{n×r} then Δ(A) = U′V^{′T} for U′ = US, V′ = V S^{−T}, where S ∈ R^{r×r} is any nonsingular matrix. Therefore another displacement representation can be derived by the latter decomposition.

If the product of the norms of the columns of U and the corresponding columns of V is much larger than the norm of A, then computing the matrix vector product Ax by means of the displacement representation formula (2.2) may lead to large numerical cancellation errors.

This drawback can be overcome by using the concept of orthogonal displacement representation.

If Δ(A) = UΣV^{T}, Σ = diag(σ_{1}, …, σ_{r}), U, V ∈ R^{n×r}, is the singular values decomposition of Δ(A), then the quadruple (a, U, σ, V) defines a set of generators, where σ = (σ_{1}, …, σ_{r}), such that A = L(a)+∑_{i=1}^{r} σ_{i}L(u_{i})L(Zv_{i})^{T}, where u_{i} = Ue_{i}, v_{i} = V e_{i}, i = 1, …, r, and e_{i} denotes the i-th column of the n × n identity matrix. We refer to this particular set of generators as to the orthogonal generators for the displacement representation of A.

Due to the uniqueness of the SVD of a matrix the set of orthogonal generators of a matrix is unique.

Observe that for the Euclidean norm ||·|| we have ||Δ(A)|| = ||AZ −ZA|| ≤ 2||A||, and for the orthogonal set of generators it holds ||U|| = ||V || = 1, σ_{1} = ||Δ(A)||, whence σ_{1} ≤ 2||A||. This fact overcomes the possibility of numerical cancellation.

If A is such that A^{T} = JAJ, like in the case of Toeplitz matrices, where J is the permutation matrix having 1 in the anti-diagonal, there exists an additional structure of Δ(A). In fact it can be easily proved that Δ(A)J is skew-symmetric.

For skew-symmetric matrices we have the following technical result:

Lemma 2.1. Let A ∈ R^{n×n} be such that A = −A^{T} and rank(A) = r. Then r = 2h is even and there exist W ∈ R^{n×r}, D ∈ R^{h×h} such that W^{T} W = I_{r}, D = diag(d_{1}, …, d_{h}) and

A = W (D ⊗ S)W^{T} , S = [ \begin{matrix}0 & 1 \\ −1 & 0 \end{matrix} ] ,

moreover ±id_{j} are the nonnull eigenvalues of A, where i is the complex unit such that i^{2} = −1, and |d_{j}| = σ_{2j−1} = σ_{2j}, j = 1, …, h, where σ_{1}, …, σ_{r} are the singular values of A.

Proof. From the Schur normal form deduce that a skew symmetric matrix A is such that A = QD̂Q^{H}, for Q^{H}Q = I, D̂ = diag(d_{j}) where d_{2j−1} = θ_{j}i, d_{2j} = −θ_{j}i, θ_{j} ∈ R, |θ_{j}| = σ_{2j−1} = σ_{2j}, j = 1, …, h, d_{j} = 0, j = 2h + 1, …, n, and A^{H} denotes the transpose conjugate of A. Therefore, only an even number of eigenvalues can be nonzero. Denote q_{j} = Qe_{j}; since A is real then q_{2j} = q̄_{2j−1}; apply the equation

[ \begin{matrix}θi & 0 \\ 0 & −θi \end{matrix} ] = \frac{θ}{2} [ \begin{matrix}1 & −i \\ 1 & i \end{matrix} ] [ \begin{matrix}0 & 1 \\ −1 & 0 \end{matrix} ] [ \begin{matrix}1 & 1 \\ i & −i \end{matrix} ]

and conclude that the claim is true where the columns of W are given by (q_{2j−1} + q̄_{2j−1})\sqrt{2} and (q_{2j−1} − q̄_{2j−1})i\sqrt{2}, that are orthonormal.

As a consequence, we have the following fundamental property:

Theorem 2.2. [Displacement representation of a persymmetric matrix]. If the persymmetric matrix A has displacement rank r with respect to Δ, then r = 2h and there exist matrices W ∈ R^{n×r}, D = diag(d_{1}, …, d_{h}) such that W^{T} W = I_{r}, |d_{1}| ≥ |d_{2}| ≥ … ≥ |d_{h}|, |d_{j}| = σ_{2j−1} = σ_{2j}, j = 1, …, h, where σ_{1}, …, σ_{r} are the singular values of Δ(A), and

Δ(A) = W (D ⊗ S)W^{T} J, S = [ \begin{matrix}0 & 1 \\ −1 & 0 \end{matrix} ] .

Therefore,

A = L(a) + ∑_{i=1}^{h} d_{i}(L(w_{2i−1})L^{T} (ZJw_{2i}) − L(w_{2i})L^{T} (ZJw_{2i−1})) (2.5)

where a = Ae_{1}, w_{i} = W e_{i}, i = 1, …, r. The orthogonal generator for A is completely defined by the triple (a, W, d) where d = (d_{i})_{i=1,...,h}.

This fact can be exploited for reducing the computational cost of multiplying the matrix A and a vector. Indeed, it is easy to prove that, given the FFT associated with x and y, for the computation of the products L(w)x and L(ZJw)y only the FFT associated with w is sufficient.

2.3. Approximate displacement rank. In certain situations it may be useful to measure the distance of a matrix A from the set of triangular Toeplitz matrices in an approximate way. In fact, it may happen that a matrix A is “very far” from the set of triangular Toeplitz matrices since its displacement rank is large (say, n), but its Euclidean distance is small, say, of the order of the machine precision. This occurs, for instance, for the matrix

A = T + E

where T is triangular Toeplitz and for the Euclidean norm ||E|| it holds ||E|| < ϵ. In this case representing A in terms of its “full” displacement representation is not suitable and is more convenient to replace A with T which is a very accurate approximation (if ϵ is small enough) and can be efficiently represented in terms of its displacement representation.

Definition 2.3. For a given ϵ > 0 define the ϵ-displacement rank of a matrix A as Drk_{ϵ}(A) = min_{||E||≤ϵ} rank(Δ(A) + E), and the relative ϵ-displacement rank as Rdrk_{ϵ}(A) = min_{||E||≤ϵ||A||} rank(Δ(A) + E).

If r = Drk_{ϵ}(A) (r = Rdrk_{ϵ}(A)), and if Δ(A) = ∑_{i=1}^{r} u_{i}v_{i}^{T} + E, with ||E|| ≤ ϵ (||E|| ≤ ϵ||A||), we say that the matrix

L(Ae_{1}) + ∑_{i=1}^{r} L(u_{i})L^{T} (Zv_{i}) (2.6)

is an ϵ-displacement (relative ϵ-displacement) representation of A.

The following properties hold.

Theorem 2.4. Let Δ(A) = UΣV^{T} be the singular value decomposition of Δ(A), where σ_{1} ≥ … ≥ σ_{k} > 0 are the singular values of Δ(A). Let σ_{1} > ϵ ≥ σ_{k}; then Drk_{ϵ}(A) = r if and only if σ_{r} > ϵ ≥ σ_{r+1}. Similarly, let σ_{1} > ϵ||A|| ≥ σ_{k}; then Rdrk_{ϵ}(A) = r if and only if σ_{r} > ϵ||A|| ≥ σ_{r+1}.

Proof. We prove that σ_{r+1} ≤ ϵ if and only if Drk_{ϵ}(A) ≤ r. This is sufficient to prove the theorem. Let u_{i}, v_{i}, i = 1, …, k, be the i-th column of U and V, respectively, that is

Δ(A) = ∑_{i=1}^{r} u_{i}σ_{i}v_{i}^{T} + ∑_{i=r+1}^{k} u_{i}σ_{i}v_{i}^{T} . (2.7)

Define E = ∑_{i=r+1}^{k} u_{i}σ_{i}v_{i}^{T}, then ||E|| = σ_{r+1}. Thus, if σ_{r+1} ≤ ϵ, we have Drk_{ϵ}(A) ≤ Drk(Δ(A) − E) = r. Vice-versa, if Drk_{ϵ}(A) = r then there exists a matrix E, with ||E|| ≤ ϵ, such that rank(Δ(A) + E) = r. Define θ_{r} = min_{{X: rank(X)=r}} ||Δ(A) − X||; since θ_{r} = σ_{r+1} and σ_{r+1} = ||Δ(A)−∑_{i=1}^{r} u_{i}σ_{i}v_{i}^{T} ||, it follows that σ_{r+1} = θ_{r} ≤ ||E||.

Theorem 2.5. Let σ_{1}, …, σ_{k} be the singular values of Δ(A), σ_{1} > ϵ ≥ σ_{k}, r = Drk_{ϵ}(A) ≤ Drk(A) = k, and let A_{ϵ} be an ϵ-displacement representation of A. Then it holds

||A − A_{ϵ}|| ≤ n ∑_{i=r+1}^{k} σ_{i} ≤ n(k − r)ϵ.

Similarly, if r = Rdrk_{ϵ}(A) and if A_{ϵ} is a relative ϵ-displacement representation of A, then it holds

\frac{||A_{ϵ} − A||}{||A||} ≤ n ∑_{i=r+1}^{k} σ_{i} ≤ n(k − r)ϵ.

Proof. Let Δ(A) = UΣV^{T} be the singular values decomposition of Δ(A). From Theorem 2.4 and (2.7) it holds A − A_{ϵ} = ∑_{i=r+1}^{k} σ_{i}L(u_{i})L(Zv_{i})^{T}, where ||u_{i}|| = ||v_{i}|| = 1. Since for any vector u it holds ||L(u)|| = ||L(u)^{T} || and ||L(u)|| ≤ ∑_{i=1}^{n} |u_{i}|· ||Z^{i}|| = ∑_{i=1}^{n} |u_{i}| ≤ \sqrt{n}||u||, we have ||A−A_{ϵ}|| ≤ n ∑_{i=r+1}^{k} σ_{i}. The proof is completed since |σ_{i}| ≤ ϵ for r = Drk_{ϵ}(A) and |σ_{i}| ≤ ϵ||A|| for r = Rdrk_{ϵ}(A), i = r + 1, …, k.

Conversely, we have the following

Theorem 2.6. If A has displacement rank k and B is an approximation of A such that ||A − B|| ≤ ϵ, then Drk_{2ϵ}(B) ≤ k. Similarly, if ||A − B|| ≤ ϵ||A|| then Rdrk_{2ϵ}(B) ≤ k.

Proof. It holds Δ(B) = Δ(A) + Δ(B − A), moreover, ||Δ(B − A)|| ≤ ||Z(B − A)|| + ||(B − A)Z|| ≤ 2||B − A|| ≤ 2ϵ since ||Z|| = 1.

Theorem 2.7. The following relations extend to the approximate displacement rank the analogous inequalities valid for displacement rank.

Drk_{δ}(A + B) ≤ Drk_{ϵ}(A) + Drk_{ϵ}(B),  δ = 2ϵ
Drk_{δ}(AB) ≤ Drk_{ϵ}(A) + Drk_{ϵ}(B),  δ = ϵ(||A|| + ||B||)
Drk_{δ}(A^{−1}) ≤ Drk_{ϵ}(A),  δ = ||A^{−1}||^{2}ϵ
Rdrk_{δ}(A + B) ≤ Rdrk_{ϵ}(A) + Rdrk_{ϵ}(B), δ = ϵ\frac{||A||+||B||}{||A+B||}
Rdrk_{δ}(AB) ≤ Rdrk_{ϵ}(A) + Rdrk_{ϵ}(B),  δ = 2ϵ\frac{||A||·||B||}{||AB||}
Rdrk_{δ}(A^{−1}) ≤ Rdrk_{ϵ}(A),  δ = ϵCond(A),

(2.8)

where Cond(A) = ||A|| · ||A^{−1}|| is the spectral condition number of A.

Proof. Let us prove the second inequality. Assuming Δ(A) = S_{A} + E_{A}, Δ(B) = S_{B} + E_{B}, where ||E_{A}|| ≤ ϵ, ||E_{B}|| ≤ ϵ, from the relation Δ(AB) = AΔ(B) + Δ(A)B, we have

Δ(AB) = AS_{B} + S_{A}B + E_{AB}, E_{AB} = AE_{B} + E_{A}B.

Thus ||E_{AB}|| ≤ ϵ(||A|| + ||B||) and rank(AS_{B} + S_{A}B) ≤ rank(S_{A}) + rank(S_{B}). Similarly we do for the remaining inequalities.

2.4. Cyclic reduction for block tridiagonal block Toeplitz systems. Cyclic reduction is a well known technique originally introduced for the solution of certain systems arising from the discretization of elliptic problems [18], and recently adjusted for the solution of Markov chains [24, 5, 6, 7]. Numerical properties of this method are described in [31, 1]. Here we recall the main properties.

Consider the system

Qx = b, (2.9)

where Q is a finite block tridiagonal block Toeplitz matrix, of block size 2^{q}, say, obtained by truncating (1.1), where Â = A. Partition the vectors x and b into blocks x_{i}, b_{i} of dimension m, respectively, that is, x = (x_{i})_{i=1,...,2^{q}}, b = (b_{i})_{i=1,...,2^{q}}. By performing an even-odd permutation of the block rows and block columns in (2.9) we find that

[ \begin{matrix}V_{0} & U_{0} \\ L_{0} & W_{0} \end{matrix} ] [ \begin{matrix}x & {(0) \atop +} \\ x & {(0) \atop −} \end{matrix} ] = [ \begin{matrix}b & {(0) \atop +} \\ b & {(0) \atop −} \end{matrix} ] (2.10)

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

\begin{matrix}V_{0} = W_{0} = [ \begin{matrix}A & & & & ◯ \\ & . & . & . \\ ◯ & & & & A \end{matrix} ] , L_{0} = [ \begin{matrix}C & & & ◯ \\ B & C \\ & ... & ... \\ ◯ & & B & C \end{matrix} ] , \\ U_{0} = [ \begin{matrix}B & C & & ◯ \\ & B & ... \\ & & ... & C \\ ◯ & & & B \end{matrix} ] . \end{matrix}

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

\{ \begin{matrix}(W_{0} − L_{0}V_{0}^{−1}U_{0})x_{−}^{(0)} = b_{−}^{(0)} − L_{0}V_{0}^{−1}b_{+}^{(0)} \\ x_{+}^{(0)} = V_{0}^{−1}(b_{+}^{(0)} − U_{0}x_{−}^{(0)}). \end{matrix} (2.11)

Denoting

Q^{(1)} = W_{0} − L_{0}V_{0}^{−1}U_{0} (2.12)

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

Q^{(1)}x^{(1)} = b^{(1)}. (2.13)

Once x^{(1)} has been computed, the solution of the original system (2.9) can be recovered by means of back substitution through (2.11).

It is easy to verify that Q^{(1)} is a block tridiagonal matrix which, except for the north western corner block, has the block Toeplitz structure. If we apply to system (2.13) the same block even-odd permutation, followed by one step of Gaussian elimination (i.e., one step of cyclic reduction), we obtain a new system, whose matrix Q^{(2)} is a block tridiagonal matrix which, except for the north western corner block, has the block Toeplitz structure. In this way we obtain a sequence of systems {Q^{(j)}x^{(j)} = b^{(j)}}, of block size 2^{q−j}, j = 0, …, q, where Q^{(0)} = Q, and

Q^{(j)} = [ \begin{matrix}Â_{j} & C_{j} & & ◯ \\ B_{j} & A_{j} & C_{j} \\ & B_{j} & A_{j} & ... \\ ◯ & & ... & ... \end{matrix} ] .

The block entries A_{j+1}, B_{j+1}, C_{j+1}, Â_{j+1}, obtained at step j + 1, that define the matrix Q^{(j+1)}, are related to the block entries A_{j}, B_{j}, C_{j}, Â_{j}, obtained at step j, by means of the following simple relations:

\{ \begin{matrix}B_{j+1} = −B_{j}A_{j}^{−1}B_{j} \\ A_{j+1} = A_{j} − B_{j}A_{j}^{−1}C_{j} − C_{j}A_{j}^{−1}B_{j} \\ C_{j+1} = −C_{j}A_{j}^{−1}C_{j} \\ Â_{j+1} = Â_{j} − C_{j}A_{j}^{−1}B_{j}, \end{matrix} (2.14)

for A_{0} = Â_{0} = A, B_{0} = B, C_{0} = C. It is useful to express these relations in functional form; for this purpose let us associate with Q^{(j)} the following matrix polynomial

ϕ_{j}(z) = C_{j} + zA_{j} + z^{2}B_{j}, j = 0, 1, 2, …. (2.15)

From (2.12) we may easily deduce that

\{ \begin{matrix}ϕ_{j+1}(z) = zA_{j} − (C_{j} + zB_{j})A_{j}^{−1}(C_{j} + zB_{j}) \\ Â_{j+1} = Â_{j} − C_{j}A_{j}^{−1}B_{j} \end{matrix}

and that

ϕ_{j+1}(z^{2}) = (\frac{ϕ_{j}(z)^{−1} − ϕ_{j}(−z)^{−1}}{2z} )^{−1} = −ϕ_{j}(z)A_{j}^{−1}ϕ_{j}(−z).

Thus, if we introduce the formal matrix power series ψ_{j}(z) = ϕ_{j}(z)^{−1}, we have

ψ_{j+1}(z^{2}) = \frac{ψ_{j}(z) − ψ_{j}(−z)}{2z}, (2.16)

that is, the matrix coefficients of ψ_{j+1}(z) are the odd matrix coefficients of ψ_{j}(z). In other words, the matrix power series ψ_{j}(z) are determined by matrix coefficients that are suitable coefficients of the matrix power series ψ(z) = ψ_{0}(z).

Relation (2.16) is useful to prove structural properties of the block entries obtained at each step of cyclic reduction. Indeed, from (2.3) it holds that for any z, ψ_{j}(z) and ϕ_{j}(z) have the same displacement rank. Unfortunately this does not imply that Drk(ψ_{j}(z)) = Drk(ψ(z)) and thus that Drk(ϕ_{j}(z)) = Drk(ψ(z)) for any j. However, under suitable conditions on ϕ(z) the former property is satisfied. This occurs for instance if ϕ(z) belongs to the same matrix algebra for any z, or even under weaker conditions. Examples related to banded Toeplitz matrices are discussed in the next section.

It is important to point out that if the matrix polynomials ϕ_{j}(z) has a small displacement rank for any j and z, then the process of cyclic reduction can be implemented at a very low cost with the help of FFT, and with a low memory storage.

Cyclic reduction has very nice convergence properties that allow one to apply this technique even to (bi-)infinite systems. In fact, it can be proved that under suitable conditions lim_{i} B_{i} = lim_{i} C_{i} = 0. More precisely we have the following result [11].

Theorem 2.8. Assume that A_{j} is nonsingular for j ≥ 0. Let ξ_{i}, i = 1, …, 2m, be the zeros of ϕ(z) = C + zA + z^{2}B and suppose that |ξ_{1}| ≤ … ≤ |ξ_{m}| < 1 < |ξ_{m+1}| ≤ … ≤ |ξ_{2m}|, where we assume that ϕ(z) has zeros equal to infinity, if its degree is less than 2m. Denote u_{1}, …, u_{2m} a set of Jordan chains associated with ξ_{1}, …, ξ_{2m}, respectively, for the matrix polynomial ϕ(z). If the vectors u_{1}, …, u_{m} are linearly independent and if the vectors u_{m+1}, …, u_{2m} are linearly independent then the sequences ||A_{j}||, ||A_{j}^{−1}|| are bounded and for any ϵ > 0 there exists an operator norm || · || such that ||C_{j}|| = O((|ξ_{m}| + ϵ)^{2^{j}}), ||B_{j}|| = O((1/|ξ_{m+1}| + ϵ)^{2^{j}}).

Proof. (Outline.) For the Jordan chains properties [23, 15] there exist m × m matrices H and K that solve the matrix equations C + HA + H^{2}B = 0 and K^{2}C + KA + B = 0 whose eigenvalues, counting multiplicities, are ξ_{1}, …, ξ_{m} and 1/ξ_{m+1}, …, 1/ξ_{2m}, respectively. It can be proved that [6] at each step j of cyclic reduction, it holds C_{j}A_{j}^{−1}+H^{2^{j}} +H^{2·2^{j}}B_{j}A_{j}^{−1} = 0, K^{2·2^{j}}C_{j}A_{j}^{−1}+K^{2^{j}} +B_{j}A_{j}^{−1} = 0. Recovering B_{j}A_{j}^{−1} from the latter equation and replacing it into the former one, yields C_{j}A_{j}^{−1} = H^{2^{j}}R_{j}, where R_{j} is a matrix with bounded norm. Since ρ(H) = |ξ_{m}| < 1, we deduce that for any ϵ there exists an operator norm || · || such that ||H|| ≤ |ξ_{m}| + ϵ, that is, ||C_{j}A_{j}^{−1}|| = O((|ξ_{m}| + ϵ)^{2^{j}}). Similarly we may do with B_{j}A_{j}^{−1}. Hence, the sequences ||A_{j}||, ||A_{j}^{−1}|| are bounded [31] and therefore ||C_{j}|| = O((|ξ_{m}| + ϵ)^{2^{j}}), ||B_{j}|| = O((1/|ξ_{m+1}| + ϵ)^{2^{j}}).

If the blocks A, B and C satisfy additional properties we can prove stronger convergence results.

Theorem 2.9. Assume that A = A^{T}, A is positive definite and one of the following two conditions holds:

1. B = C^{T}, γ(z) = z̄C + A + zB is positive definite for |z| = 1, z ∈ C;

2. B = −C^{T}.

Then:

– for j ≥ 1 B_{j} = C_{j}^{T} and A_{j}, Â_{j} are positive definite;

– ||A_{j}||, ||A_{j}^{−1}||, ||Â_{j}|| are bounded from above by a constant;

– for any ϵ > 0 there exists a norm || · || such that ||C_{j}|| = O((σ + ϵ)^{2^{j}}), where σ = max{|ξ| : det γ(ξ) = 0, |ξ| < 1}.

Moreover, in the hypothesis 1 also ||Â_{j}^{−1}|| is bounded.

Proof. Assume that condition 1 holds. Consider the block tridiagonal matrix T_{n} of block size n obtained by truncating (1.1) at block size n with Â = A. Since the eigenvalues of T_{n} belong to the interval [λ_{1}, λ_{2}], where 0 < λ_{1} = min_{|z|=1} λ(γ(z)), λ_{2} = max_{|z|=1} λ(γ(z)) [30], and λ(.) denotes the set of eigenvalues, then T_{n}, as well as any Schur complement in T_{n}, is positive definite. Therefore, A_{j}, Â_{j}, being suitable Schur complements in T_{n} for suitable values of n [10, 4], have eigenvalues in [λ_{1}, λ_{2}] so that they are positive definite and ||A_{j}||, ||A_{j}^{−1}||, ||Â_{j}||, ||Â_{j}^{−1}|| are bounded. Under hypothesis 1, there exist solutions X, Y to the equations B^{T} Y^{2} − AY + B = 0, BX^{2}−AX +B^{T} = 0 such that for the spectral radius ρ(·) it holds ρ(X) = ρ(Y) = σ < 1 [25]. By applying cyclic reduction we obtain that [6, 7] B_{j}^{T} Y^{2·2^{j}} − A_{j}Y^{2^{j}} + B_{j} = 0, B_{j}X^{2·2^{j}} − A_{j}X^{2^{j}} + B_{j}^{T} = 0. Recovering B_{j}^{T} from the latter and substituting it into the former, we arrive at

B_{j}(I − X^{2·2^{j}}Y^{2·2^{j}}) = A_{j}(I − X^{2^{j}}Y^{2^{j}})Y^{2^{j}}.

Since ρ(X) = ρ(Y) < 1 and ||A_{j}|| is bounded, for any ϵ > 0 there exists a norm || · || such that ||B_{j}|| = O((σ + ϵ)^{2^{j}}).

In the assumption 2, since A is nonsingular we may perform one step of cyclic reduction. The matrices A_{1}, Â_{1}, B_{1} C_{1} that we obtain in this way are such that B_{1} = C_{1}^{T}, A_{1}, Â_{1} are symmetric and positive definite. Moreover, the function γ_{1}(z) = zB_{1} + A_{1} + z̄B_{1}^{T} is such that γ_{1}(z) = A + (B − z̄B^{T})A^{−1}(B^{T} − zB), that is positive definite for |z| = 1. The thesis follows from part 1.

3. Solving specific structured problems. In this section we analyze the cyclic reduction algorithm, applied to certain block tridiagonal block Toeplitz matrices, having particular block entries.

3.1. The case of banded Toeplitz matrices: application to polynomial factorization. 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}. Let Q_{n} = (q_{i,j}) the matrix of size n having entries q_{i,j} = a_{i−j}. Q_{n} is a banded Toeplitz matrix, with bandwidth m = max{m_{1}, −m_{2}}. Assuming that n = m2^{q}, the matrix Q_{n} can be partitioned into m × m blocks yielding a block tridiagonal Toeplitz matrix with the structure (1.1), where

\begin{matrix}Â = A = [ \begin{matrix}a_{0} & a_{−1} & … & a_{−m+1} \\ a_{1} & a_{0} & ... & ⋮ \\ ⋮ & ... & ... & a_{−1} \\ a_{m−1} & … & a_{1} & a_{0} \end{matrix} ] , \\ C = [ \begin{matrix}a_{−m} & 0 & … & 0 \\ a_{−m+1} & a_{−m} & ... & ⋮ \\ ⋮ & ... & ... & 0 \\ a_{−1} & … & a_{−m+1} & a_{−m} \end{matrix} ] , B = [ \begin{matrix}a_{m} & a_{m−1} & … & a_{1} \\ 0 & a_{m} & ... & ⋮ \\ ⋮ & ... & ... & a_{m−1} \\ 0 & … & 0 & a_{m} \end{matrix} ] . \end{matrix} (3.1)

In order to solve Q_{n}x = b, we apply the cyclic reduction algorithm to the block tridiagonal system, as described in section 2.4. The Toeplitz structure of the blocks A, B, C is generally lost by the blocks A_{j}, B_{j}, C_{j}, Â_{j}, obtained in the cyclic reduction process. This fact might lead to think that the cyclic reduction algorithm does not fully exploit the Toeplitz structure of the matrix Q. However it is possible to prove that the displacement rank of the blocks A_{j}, B_{j}, C_{j}, Â_{j} is at most 4. This fact has been proved in [10]. Here we provide a simpler proof of this property that gives more insight to possible extensions and generalizations of it.

Indeed, it is easy to verify that for any z the matrix ϕ(z) = C + zA + z^{2}B is a z-circulant matrix, i.e., belongs to the algebra generated by

C_{z} = [ \begin{matrix}0 & & & z \\ 1 & 0 \\ & ... & ... \\ & & 1 & 0 \end{matrix} ] .

Since this set is closed under multiplication and addition, also the inverse ψ(z) = ϕ(z)^{−1} is a z-circulant matrix, where ϕ(z) is nonsingular. Moreover, since z-circulant matrices are Toeplitz, it follows that ψ(z) is Toeplitz for any z. This implies that the matrices ψ_{j}(z) are Toeplitz for any z for which they exist, due to (2.16). Therefore ϕ_{j}(z) = ψ_{j}(z)^{−1}, has displacement rank at most 2 for any j and for any z for which ϕ_{j}(z) exists, due to (2.3). In particular we deduce that C_{j} = ϕ_{j}(0) has displacement rank at most 2, A_{j} = (ϕ_{j}(1) + ϕ_{j}(−1))/2 has displacement rank at most 4, and B_{j} = z^{2}ϕ_{j}(z^{−1})|_{z=0} has displacement rank at most 2.

In [10] we provide efficient formulae for the computation of the displacement generators of A_{j}, B_{j}, C_{j}, Â_{j}, that lead to O(m log^{2} m) ops the cost for performing one step of cyclic reduction.

For the polynomial factorization problem described in section 1.2, we have to compute the vector w^{(n)}, made up by the first m entries of the first column of the inverse of the matrix Q_{n}, for a large enough value of n, where Q_{n} is the n × n banded Toeplitz matrix Q_{n} = (q_{i,j}), such that q_{i,j} = p_{i−j+m} if 0 ≤ i − j + m ≤ 2m, q_{i,j} = 0, otherwise. Thus, if n = m2^{q}, we may reblock the matrix Q_{n} into a block tridiagonal block Toeplitz matrix, with blocks having the structure (3.1), where a_{i} = p_{i+m}, i = −m, …, m.

In [4] we show that, if we compute the sequence (2.14), j = 1, …, q, then it holds Â_{j}^{−1}e_{1} = w^{(m2^{j})}. Thus, formulae (2.14) provide the vector w^{(n)}, for n = m2^{j}, j = 1, 2, 3, …, and, for the displacement structure of the matrices Â_{j}, A_{j}, B_{j}, C_{j}, the cost of switching from n = m2^{j} to n = m2^{j+1} is O(m log^{2} m) arithmetic operations. Moreover, applying cyclic reduction algorithm leads to a quadratically convergent algorithm, since f = −Uw^{(m2^{j})} + O(σ^{m2^{j}}), j = 1, 2, …, where σ = |z_{m}/z_{m+1}|.

We may increase the numerical performances of the algorithm by means of the orthogonal displacement representations. In fact it easy to check that A_{j}, B_{j} and C_{j} are persymmetric for any j, thus by using (2.14) and Theorem 2.2 it is easy to compute the orthogonal displacement generators of A_{j+1}, B_{j+1}, C_{j+1} and Â_{j+1} given the orthogonal displacement generators of A_{j}, B_{j}, C_{j}, and Â_{j}.

For a persymmetric matrix X let us denote ODG(X) = (x, W_{X}, d_{X}), the orthogonal displacement generator of X as given in Theorem 2.2, and let us define ODG(X) = (x, U_{X}, σ_{X}, V_{X}), for a general matrix X, where U_{X}Σ_{X}V_{X}^{T} is the SVD of Δ(X), Σ_{X} = diag(σ_{X}) and x is the first column of X. Moreover set D_{X} =

diag(d_{X}) ⊗ S where S = [ \begin{matrix}0 & 1 \\ −1 & 0 \end{matrix} ].

Theorem 3.1. [Displacement generator in cyclic reduction]. Let us denote A′, B′, C′ and Â′ the blocks obtained by applying one step of CR to the blocks A, B, C and Â, where A, B, C are persymmetric. Then, we have:

Δ(A′) = F_{A} [ \begin{matrix}D_{A} & O & O \\ O & O & D′ \\ O & D′ & O \end{matrix} ] F_{A}^{T} J, D′ = [ \begin{matrix}−D_{B} & O & O \\ O & D_{A} & O \\ O & O & −D_{C} \end{matrix} ]

F_{A} = [W_{A}, W_{B}, BA^{−1}W_{A}, BA^{−1}W_{C}, CA^{−1}W_{B}, CA^{−1}W_{A}, W_{C}]

Δ(B′) = F_{B} [ \begin{matrix}D_{A} & O & O \\ O & O & −D_{B} \\ O & −D_{B} & O \end{matrix} ] F_{B}^{T} J, F_{B} = [BA^{−1}W_{A}, W_{B}, BA^{−1}W_{B}]

Δ(C′) = F_{C} [ \begin{matrix}D_{A} & O & O \\ O & O & −D_{C} \\ O & −D_{C} & O \end{matrix} ] F_{C}^{T} J, F_{C} = [CA^{−1}W_{A}, W_{C}, CA^{−1}W_{C}]

Δ(Â′) = F_{A} diag(Σ_{A}, −D_{C}, D_{A}, −D_{B})GT ; _{A}

F_{A} = [U_{A}, W_{C}, CA^{−1}W_{A}, CA^{−1}W_{B}] , G_{A} = [V_{A}, JBA^{−1}W_{C}, JBA^{−1}W_{A}, JW_{B}] .

The above theorem provides a displacement generator of A, B, C, Â of the form UΣV^{T} where the number of columns of U and V is generally much larger than the displacement rank. To compress this redundant displacement representation we may apply the following

Algorithm 3.2. Computing the orthogonal generator for a persymmetric matrix.

Input: Matrices F ∈ R^{m×k}, M ∈ R^{k×k} such that Δ(A) = F MF^{T} J, M = −M^{T};

k ≥ h, where rank Δ(A) = h. Output: The orthogonal generator of A in the form (a, W, d), where W ∈ R^{h×h}. Computation:

Stage 1: Compute the QR factorization U = QR.

Stage 2: Compute H = RMR^{T}.

Stage 3: Compute the SVD H = UΣV^{T}, U, V ∈ R^{k×h}, Σ ∈ R^{h×h}.

Stage 4: Compute the sign γ_{j} = sign(u_{i,2j−1}v_{i,2j}), j = 1, …, h, where i is

such that the product is nonzero. Stage 4: Set W = QU, d_{i} = γ_{i}σ_{2i−1}, i = 1, …, h, d = (d_{i}).

For a general matrix A the above algorithm must be modified in the following way

Algorithm 3.3. Computing the orthogonal generator for a general matrix.

Input: Matrices F_{1}, F_{2} ∈ R^{m×k}, M ∈ R^{k×k} such that Δ(A) = F_{1}MF_{2}^{T}, k ≥ h,

h = rank Δ(A). Output: The orthogonal generator of A in the form (a, U, σ, V), U, V ∈ R^{m×h}. Computation:

Stage 1: Compute the QR factorizations F_{i} = Q_{i}R_{i}, i = 1, 2.

Stage 2: Compute H = R_{1}MR_{2}^{T}.
Stage 3: Compute the SVD H = G_{1}ΣG_{2}^{T}, G_{1}, G_{2} ∈ R^{k×h}.
Stage 4: Set U = Q_{1}G_{1}, V = Q_{2}G_{2}, σ = (σ_{i})_{i=1,h}.

In the cyclic reduction algorithm the inversion of the Toeplitz-like matrix A_{i} must be performed. This can be obtained by means of superfast algorithms in O(m log^{2} m) ops or by means of the preconditioned conjugate gradient techniques. In both cases the numerical stability of performing the matrix-vector product is crucial for the reliability of the overall computation. Using orthogonal displacement representations may strongly improve numerical stability. An alternative way for Toeplitz-like matrix inversion is described in the last section of this paper and relies on Newton’s iteration and on the approximate displacement rank. Even in this case, orthogonal displacement representations are fundamental for the stability of this computation. Using orthogonal displacement representations does not substantially increase the cost of the overall computation.

3.2. The case of PH/PH/1 queueing problems. For the matrix of (1.1) with blocks (1.2), we may prove that [3] at each step j ≥ 1 of cyclic reduction it holds

A_{j} = A − (I_{m_{1}} ⊗ S^{0})σ_{j}(α ⊗ I_{m_{2}}) − (T^{0} ⊗ I_{m_{2}})θ_{j}(I_{m_{1}} ⊗ β)
B_{j} = −(I_{m_{1}} ⊗ S^{0})δ_{j}(I_{m_{1}} ⊗ β)
C_{j} = −(T^{0} ⊗ I_{m_{2}})γ_{j}(α ⊗ I_{m_{2}})
Â_{j} = A − (T^{0} ⊗ I_{m_{2}})θ_{j}(I_{m_{1}} ⊗ β),

where δ_{j}, γ_{j}, σ_{j}, θ_{j} are m_{1} × m_{1}, m_{2} × m_{2}, m_{1} × m_{2}, m_{2} × m_{1} matrices, respectively. The matrices δ_{j}, γ_{j}, σ_{j}, θ_{j}, for j ≥ 2, are defined by the recursions

δ_{j+1} = δ_{j} (δ_{1} + [ δ_{1} σ_{1} ] D_{j}V_{j}^{−1} [ \begin{matrix}θ_{1} \\ δ_{1} \end{matrix} ]) δ_{j},

γ_{j+1} = γ_{j} (γ_{1} + [ θ_{1} γ_{1} ] D_{j}V_{j}^{−1} [ \begin{matrix}γ_{1} \\ σ_{1} \end{matrix} ]) γ_{j},

σ_{j+1} = σ_{j} + δ_{j} (σ_{1} + [ δ_{1} σ_{1} ] D_{j}V_{j}^{−1} [ \begin{matrix}γ_{1} \\ σ_{1} \end{matrix} ]) γ_{j},

θ_{j+1} = θ_{j} + γ_{j} (θ_{1} + [ θ_{1} γ_{1} ] D_{j}V_{j}^{−1} [ \begin{matrix}θ_{1} \\ δ_{1} \end{matrix} ]) δ_{j},

for j = 1, 2, …, where

δ_{1} = (I_{m_{1}} ⊗ β)A^{−1}(I_{m_{1}} ⊗ S^{0}),
γ_{1} = (α ⊗ I_{m_{2}})A^{−1}(T^{0} ⊗ I_{m_{2}}),
σ_{1} = (I_{m_{1}} ⊗ β)A^{−1}(T^{0} ⊗ I_{m_{2}}),
θ_{1} = (α ⊗ I_{m_{2}})A^{−1}(I_{m_{1}} ⊗ S^{0}),

and

D_{j} = [ \begin{matrix}σ_{j} & 0 \\ 0 & θ_{j} \end{matrix} ] , V_{j} = I − [ \begin{matrix}θ_{1} & γ_{1} \\ δ_{1} & σ_{1} \end{matrix} ] D_{j}.

In this way instead of processing matrices of size m_{1}m_{2}, we have to deal just with matrices of size m_{1} + m_{2}. The cost of one step is so reduced from O(m_{1}^{3}m_{2}^{3}) to O(m_{1}^{3} + m_{2}^{3}) ops. The algorithm has been implemented and the results of the numerical tests proved its reliability and efficiency. More details can be found in [3].

3.3. The Ramaswami reduction: Frobenius structures. One of the main problems in the solution of M/G/1 type Markov chains is the computation of the minimal nonnegative solution G of the nonlinear matrix equation G = ∑_{i=0}^{∞} D_{i}G^{i}, where D_{i}, i ≥ 0, are m × m nonnegative matrices such that ∑_{i=0}^{∞} D_{i} is stochastic. In [28], V. Ramaswami has proved that this problem is reduced to the problem of solving the infinite block tridiagonal system

Q = [ \begin{matrix}Â & C & & ◯ \\ B & A & C \\ & B & A & ... \\ ◯ & & ... & ... \end{matrix} ] [ H ; H^{2} ; H^{3} ; ⋮ ] = [ B ; 0 ; 0 ; ⋮ ] , (3.2)

with infinite blocks defined by

\begin{matrix}A = [ \begin{matrix}I − D_{1} & −D_{2} & −D_{3} & … \\ & I \\ & & I \\ . & . & . \end{matrix} ] , B = [ \begin{matrix}−D_{0} & 0 & … \\ 0 & 0 & … \\ ⋮ & ⋮ \end{matrix} ] , \\ C = [ \begin{matrix}0 & 0 & … \\ −I & 0 & ... \\ & −I & ... \\ . & . & . \end{matrix} ] . \end{matrix} (3.3)

In [12] we show that, if D_{i} = 0, for i > N, then the system (3.2) can be reduced to an infinite block tridiagonal system, with N × N blocks, obtained by taking finite N-dimensional sections of matrices (3.3). Moreover, under the notation of (2.14) if we apply cyclic reduction to the latter system, the sequence W_{j} = −(I − Â_{j})^{−1}B is such that H = W_{j} +O(σ^{2^{j}}) where 0 < σ < 1, and where A_{0}, B_{0}, C_{0} are the N-dimensional finite sections of the matrices (3.3), and H is the N-dimensional finite section of the corresponding matrix in (3.2). The efficient computation of the blocks generated by the cyclic reduction process relies on the Frobenius structure of A_{j}, B_{j}, C_{j}. Indeed, it is a simple matter to show that the matrices A_{j} and B_{j} remain structured at each step j: B_{j} is a matrix having all the entries equal to zero, except for the entries of the block (1,1); A_{j} has the structure I + u_{j}e_{1}^{T} + e_{1}v_{j}^{T}, where v_{j}, u_{j} are suitable block vectors, and e_{1} is the first block column of the identity matrix. On the other hand the matrix C_{j}, after log N steps is a full matrix, without any evident structure. In fact, also the matrix C_{j} is structured, and the structure is related to the algebra generated by a block Frobenius matrix. More precisely, if D_{0} is nonsingular, define the block Frobenius matrix

F = [ \begin{matrix}−D_{0}^{−1}(D_{1} − I) & −D_{0}^{−1}D_{2} & … & −D_{0}^{−1}D_{N} \\ I & 0 & … & 0 \\ & ... & ... & ⋮ \\ ◯ & & I & 0 \end{matrix} ] .

Then the matrix polynomial ϕ(z) = B + zA + z^{2}C can be written as

ϕ(z) = [ \begin{matrix}D_{0} & & & & & ◯ \\ & −zI \\ & & . & . & . \\ ◯ & & & & & −zI \end{matrix} ] (I − zF)

that is, ϕ(z) is a block diagonal matrix times a matrix polynomial with block coefficients defined by powers of F. A similar structure is kept at each step j; more precisely, we prove that the matrix polynomial of (2.15) ϕ_{j}(z) = C_{j} + zA_{j} + z^{2}B_{j} can be expressed as

ϕ_{j}(z) = (I − w_{j}(z)H_{j}e_{1}^{T})(zI − F^{2^{j}}) (3.4)

where H_{j} = (e_{1}^{T} F^{2^{j}−1}e_{1})^{−1}, w_{j}(z) = F^{2^{j}−1}e_{1} + ze_{1}D_{0}. That is, up to a small rank correction, ϕ_{j}(z) is a matrix polynomial with block coefficients generated by powers of F. From this property it follows that B_{j} = ϕ_{j}(0) has displacement rank at most 2. By using a continuity argument, we prove that the same bound holds also in case where D_{0} is singular.

Thus, the computation of A_{j}, B_{j}, C_{j} costs O(m^{3}N + m^{2}N log N) ops.

3.4. The case of block banded block Toeplitz matrices. Consider the block banded block Toeplitz matrix Q with banded Toeplitz blocks, described in section 1.3. By partitioning Q into (mM) × (mM) blocks we obtain a block tridiagonal matrix with blocks A, B, C. It is immediate to verify that these blocks have displacement rank at most 2M. In particular, for the matrix polynomial ϕ(z) and for the matrix function ψ(z) = ϕ(z)^{−1} it holds Drk(ϕ(z)) = Drk(ψ(z)) ≤ 2M, for any z.

The matrix polynomial ϕ(z) has no algebraic properties that could lead to a conservation of the displacement rank. Even in the case M = 1, where ϕ(z) is a tridiagonal Toeplitz matrix, the matrix polynomials ϕ_{j}(z) do not maintain a constant displacement rank.

A possible way for overcoming this drawback is to look for different displacement operators associated with different matrix algebras like trigonometric algebras. This is a research direction to be better investigated.

A different approach is based on the concept of approximate displacement rank introduced in section 2.3. Indeed, under the assumptions that guarantee the convergence to zero of the blocks C_{j}, B_{j} (compare with Theorems 2.8, 2.9), it is easy to prove that for a given ϵ, Drk_{ϵ}(A_{j}) is bounded by a constant independent of j. Similarly it holds for Rdrk_{ϵ}(B_{j}), and Rdrk_{ϵ}(C_{j}).

The following modification of cyclic reduction that generates approximations Ã_{j},

B̃_{j}, C̃_{j}, ˜ Â_{j}, of the blocks A_{j}, B_{j}, C_{j}, Â_{j}, expressed in terms of their relative ϵdisplacement generator is based on the above property.

For simplicity we just give an outline of the algorithm in its general j-th step.

Algorithm 3.4. Input: A relative ϵ-displacement generator of Ã_{j}, B̃_{j}, C̃_{j}, ˜ Â_{j}.

Output: A relative ϵ-displacement generator of Ã_{j+1}, B̃_{j+1}, C̃_{j+1}, ˜ Â_{j+1}. Computation:

1. Compute the first column b̃ = −B̃Ã^{−1}B̃e_{1}.

2. By using the relative ϵ-displacement generator of Ã_{j}, compute a relative δdisplacement generator of Ã_{j}^{−1}, for δ = ϵCond(Ã_{j}).

3. From the latter, and from the relative ϵ-displacement generator of B̃_{j} compute a relative θ-displacement generator of

H = −B̃_{j}Ã_{j}^{−1}B̃_{j},

for a suitable θ, according to formulae (2.8). Let h be the rank of such generator and S_{1} and S_{2} be the mM × h matrices whose columns define a generator, i.e., such that

Δ(−B̃_{j}Ã_{j}^{−1}B̃_{j}) = S_{1}S_{2}^{T} + E, ||E|| ≤ θ||B̃_{j}Ã_{j}^{−1}B̃_{j}||.

4. Compute the QR factorizations S_{1} = Q_{1}R_{1}, S_{2} = Q_{2}R_{2} and the h×h matrix W = R_{1}R_{2}^{T}. Compute the singular value decomposition W = UΣV^{T} of W and denote σ_{1}, …, σ_{h} the singular values of W. Let r be such that σ_{r} > ϵ||B̃_{j}Ã_{j}^{−1}B̃_{j}|| ≥ σ_{r+1} and denote with Û, V̂, the matrices made by the first r columns of U and V, respectively. Similarly denote with Σ̂ the leading r × r principal submatrix of Σ.

5. Define the matrix B̃_{j+1} as the one represented by the displacement generator (b̃, Q_{1}Û, σ̂, Q_{2}V̂), σ̂ = (σ_{1}, …, σ_{r}), i.e., B̃_{j+1} is the matrix such that Δ(B̃_{j+1}) = (Q_{1}Û)Σ̂(Q_{2}V̂)^{T}. This matrix provides a relative ϵ-displacement representation of H.

6. Perform similar computations for the matrices −C̃_{j}Ã_{j}^{−1}C̃_{j}, Ã_{j} − B̃_{j}Ã_{j}^{−1}C̃_{j} −

C̃_{j}Ã_{j}^{−1}B̃_{j} and for ˜ Â_{j} − B̃_{j}Ã_{j}^{−1}C̃_{j}.

The cost of the algorithm is O(h^{2}mM log^{2}(mM)) ops per step and strongly depends on the largest values of h encountered in the computation. The algorithm has been implemented and tested with certain PSF’s like the one defined by M = 4, α_{0,0} = γ/σ, α_{i,j} = 1/σ otherwise, with the normalizing factor σ = (M + 1)^{2} − 1 + γ. For γ > 2 the algorithm has been capable to reconstruct the image with a good accuracy even with rather large values of ϵ, say in the range [10^{−6}, 10^{−2}].

Issues like the numerical stability and the control of the level ϵ where the singular values are truncated are still to be investigated.

The following tables report the growth of the (relative) displacement rank of the blocks A_{j}, B_{j}, C_{j} and Â_{j} generated by Algorithm 3.4 when applied to the block 5-diagonal block Toeplitz matrix having Toeplitz blocks with 5 diagonals, where all the nonzero off-diagonal entries are 1 and the diagonal entry is α.

ϵ = 10^{−2}  ϵ = 10^{−4}  ϵ = 10^{−6}
m Drk_{ϵ} Rdrk_{ϵ} Drk_{ϵ} Rdrk_{ϵ} Drk_{ϵ} Rdrk_{ϵ}
64 12  12  20  22  28  34
128 12  12  20  22  28  34
256 12  12  20  22  28  34

Table 3.1 α = 8

\begin{matrix}m & Drk_{ϵ} & Rdrk_{ϵ} \\ 64 & 52 & 22 \\ 128 & 88 & 42 \\ 256 & 150 & 62 \end{matrix}

Table 3.2 α = 2, ϵ = 10^{−2}

In tables 3.1, 3.2 we report the size m of the blocks of the reblocked tridiagonal matrix and the largest values of Drk_{ϵ} and Rdrk_{ϵ} of the blocks generated by CR, for different values of ϵ. We observe that for α = 2 the hypotheses of theorem 2.8 are not satisfied and even though the Rdrk_{ϵ} depends on N, its growth is rather moderate.

3.5. Inverting Toeplitz-like matrices with Newton’s iteration. Newton’s iteration is a valid tool for solving matrix problems where the solution to be computed can be formally viewed as the zero of a known function. For instance, the formal application of this iteration to the matrix equation X^{−1} − A = 0 yields

X_{i+1} = 2X_{i} − X_{i}AX_{i}, i ≥ 0 (3.5)

that for a suitable X_{0} provides a sequence converging to A^{−1} if det A ≠ 0. This method has been analyzed in [13] and, concerning Toeplitz matrices, has been investigated in [27] where the iteration (3.5) has been modified by replacing each X_{i} with the closest matrix having displacement rank 2. Conditions on X_{0} have been given in order that the method, modified in this way, still converges.

Observe that for a nonsingular n × n Toeplitz matrix A, its inverse A^{−1} has displacement rank at most 2, whereas the matrices X_{i} generated by Newton’s iteration have displacement rank that generally grows and after few steps coincides with n. Despite the displacement rank of X_{i} is large, the Euclidean distance of X_{i} from the set of matrices with displacement rank at most 2 tends quadratically to zero. This is a crucial fact that suggests us to apply the concept of approximate displacement rank.

Here we recall the different approach presented in [9], for inverting a Toeplitzlike matrix A having displacement rank k, based on the idea of adaptivity, where each approximation X_{i} is replaced by the closest matrix having displacement rank k_{i}, and k_{i} is chosen so that the Euclidean distance of X_{i+1} from the matrices with displacement rank k is “sufficiently” small in order to guarantee the convergence.

In this way, the displacement rank is not bounded a priori but it may vary according to the difficulty of the problem and to the needs of the computation. This approach is adaptive in the sense that the cost of the algorithm is not uniform with respect to the set of input data, where the algorithm adjusts the computation to the specific difficulty of the problem. In practice, for well conditioned matrices the growth of the ϵ-displacement rank is negligible and the overall cost is O(n log n) ops. For (polynomially) ill-conditioned matrices the approximate displacement rank seems to grow as O(log n) so that the overall cost is O(n log^{2} n) ops.

We refer the reader to [9] for all the issues concerning convergence, preconditioning, regularization properties and generalizations of Newton’s iterations. Here we just want to describe how approximate displacement rank can be combined with Newton’s iteration in order to increase the performances of this method if applied to Toeplitz-like matrices.

Let us assume that Drk(A) = k and consider the following modification of (3.5)

\begin{matrix}Y_{i} = 2X_{i} − X_{i}AX_{i} \\ X_{i+1} = Trunc_{ϵ_{i}}(Y_{i}), i ≥ 0, \end{matrix} (3.6)

where Trunc_{ϵ_{i}}(Y_{i}) is the relative ϵ_{i}-displacement generator of Y_{i} given according to (2.6), and ϵ_{i} is a parameter that is dynamically computed during the iteration. The value of ϵ_{i} should be not too small to allow the displacement rank of X_{i} to reach large values, and, at the same time, should be not too large to spoil the convergence of the sequence X_{i}.

The choice of ϵ_{i} can be made according to different strategies. The simplest one is to choose a constant value of ϵ_{i}, say, the square root of the machine precision of the computation. A dynamic choice of ϵ that guarantees the convergence of the modified Newton method can be easily designed. We report below a convergence theorem, here adjusted to the operator Δ(A) = ZA − AZ and originally proved in [9] for a slightly different operator.

Theorem 3.5. Let X_{0} be such that ||X_{0}|| ≤ 1−θ, set ϵ_{i} = θ/(4n(2h_{i}+k)σ_{1}^{(i)}||A||), h_{i} = Drk(X_{i}), σ_{1}^{(i)} = ρ(Δ(X_{i})), k = Drk(A) . Then for the residuals R_{i} = I −AX_{i} of the sequence {X_{i}}_{i} generated by the modified Newton iteration (3.6) it holds ||R_{i}|| ≤ (1 − θ/2)^{2^{i}}.

For the computation of a single step of modified Newton’s iteration observe that

Δ(Y_{i}) = 2Δ(X_{i}) − Δ(X_{i})AX_{i} − X_{i}Δ(A)X_{i} − X_{i}AΔ(X_{i}) (3.7)

so that Drk(Y_{i}) ≤ k + 2Drk(X_{i}) where k = Drk(A).

Moreover, replacing for notational simplicity X_{i} with X and Y_{i} with Y, if Δ(X) = U_{X}D_{X}V_{X}^{T}, Δ(A) = U_{A}D_{A}V_{A}^{T}, where U_{X}, V_{X}, U_{A}, V_{A} ∈ R^{n×h} are orthogonal, we obtain that

\begin{matrix}Δ(Y) = U_{Y} D_{Y} V_{Y}^{T} \\ U_{Y} = [U_{X}, XU_{A}, XAU_{X}], V_{Y} = [V_{X}, X^{T} V_{A}, X^{T} A^{T} V_{X}] \\ D_{Y} = [ \begin{matrix}2D_{X} & 0 & −D_{X} \\ 0 & −D_{A} & 0 \\ −D_{X} & 0 & 0 \end{matrix} ] \end{matrix} (3.8)

Thus, an orthogonal generator of X_{i+1} = Trunc_{ϵ}(Y) can be computed according to the following

Algorithm 3.6. Computing the approximate orthogonal generator of Y = 2X − XAX.

1. Compute the first column of X_{i+1}, i.e., 2X_{i}e_{1} − X_{i}AX_{i}e_{1};
2. compute the matrices U_{Y}, V_{Y}, D_{Y} of (3.8);
3. compute the QR factorizations U_{Y} = Q_{1}R_{1}, V_{Y} = Q_{2}R_{2};

4. compute the SVD R_{1}D_{Y} R_{2}^{T} = UΣV, and r such that for the singular values of R_{1}D_{Y} R_{2}^{T} it holds σ_{r+1} ≤ ϵσ_{1} < σ_{r};

5. set U_{X} = Q_{1}Û, V_{X} = Q_{2}V̂, d_{X} = (σ_{1}, …, σ_{r}), where Û and V̂ are made up by the first r columns of U and V, respectively.

The cost of the above computation is about O(hk) FFTs. Improvements of this computation can be obtained by suitably using further structures, like the persymmetry (as in the case of Toeplitz matrices) or the symmetry of the matrices involved, that allow the reduction of the number of FFTs.

Below we resume our modification of Newton’s iteration for the inversion of a Toeplitz-like matrix assigned in terms of its orthogonal displacement generator.

Algorithm 3.7. Newton’s iteration with ϵ-displacement rank.

Input Integers n, k, the maximum number MAXIT of iterations, a residual error bound ϵ, the orthogonal generator (a, U_{A}, σ_{A}, V_{A}) of the nonsingular matrix A ∈ R^{n×n} having displacement rank k; the orthogonal generator (x, U_{X}, σ_{X}, V_{X}) of a matrix X_{0} such that AX_{0} is positive definite and ρ(AX_{0}) < 1; a strategy for the choice of the truncation value ϵ_{i}.

Output An approximation X of A^{−1} such that ||I − AX|| ≤ ϵ given in terms of its

orthogonal generator (x, U_{X}, σ_{X}, V_{X}). Computation

Stage 0 Set i = 0, X = X_{0};

Stage 1 Compute the truncation value ϵ_{i};

Stage 2 Compute the approximate orthogonal generator (y, U_{Y} , σ_{Y} , V_{Y}) of Y = 2X − XAX, by means of algorithm 3.6 with ϵ = ϵ_{i} and set (x, U_{X}, σ_{X}, V_{X}) = (y, U_{Y} , σ_{Y} , V_{Y})

Stage 3 Compute ||R|| for the residual R = I − AX; if ||R|| > ϵ and i ≤

MAXIT set i = i + 1 and continue from Stage 1. Stage 4 If i > MAXIT then output FAILURE, otherwise output the or-

thogonal generator of X.

Concerning the choice of ϵ_{i}, Theorem 3.5 provides the strategy that guarantees convergence. Different strategies based on heuristics can be applied. The starting matrix X_{0} can be chosen according to the features of A: for a general A, we may set X_{0} = A^{T} /||AA^{T} ||; for A symmetric and positive definite we may set X_{0} = I/||A||. A more convenient choice can be done in terms of the least square circulant matrix that approximates A.

3.5.1. Numerical experiments. We have applied the modified Newton iteration for inverting some test Toeplitz matrices. In Tables 3.3, 3.4, 3.5, 3.6 we report the number of steps Nsteps, the maximum displacement rank Mdrk, the sum of the displacement ranks Sdrk of each step, and the residual error err for inverting the matrix A_{n}, for each value of n, where A_{n} = (a_{|i−j|})_{i,j=1,n} is defined in the following way: Matrix 1 of Table 3.3: a_{|i−j|} = 1/(|i − j| + 1); Matrix 2 of Table 3.4: a_{0} = 4, a_{1} = 1, a_{i} = 0 otherwise; Matrix 3 of Table 3.5: a_{0} = 2.1, a_{1} = 1, a_{i} = 0 otherwise; Matrix 4 of Table 3.6: a_{0} = 2, a_{1} = 1, a_{i} = 0 otherwise. The quantity Sdrk · n log n expresses the overall arithmetic cost of the algorithm.

The matrices 1,2,3 are asymptotically well conditioned while the matrix 4 has condition number that grows as O(n^{2}).

n  Nstep Mdrk Sdrk  err
32  8  3  16  2.4e-14
64  8  3  16  5.1e-14
128  8  3  16  6.6e-14
256  8  3  16  6.8e-14
512  8  3  16  6.0e-14
1024  8  3  16  5.8e-14
2048  8  3  16  5.3e-14
4096  8  3  16  4.7e-14

Table 3.3

n  Nstep Mdrk Sdrk  err
32  8  3  15  2.0e-16
64  8  3  15  4.0e-16
128  8  3  15  2.0e-16
256  8  3  15  2.0e-16
512  8  3  15  2.0e-16
1024  8  3  15  2.1e-16
2048  8  3  15  2.0e-16
4096  8  3  15  2.2e-16

Table 3.4

n  Nstep Mdrk Sdrk  err
32  11  4  38  3.3e-14
64  12  5  41  1.7e-15
128  11  8  52  2.0e-15
256  11  7  47  2.1e-15
512  10  6  41  3.4e-13
1024  10  5  38  2.5e-13
2048  10  5  35  1.9e-13
4096  10  4  29  2.2e-13

Table 3.5

n  Nstep Mdrk Sdrk  err
32  12  5  47  3.8e-15
64  13  7  57  2.5e-15
128  14  8  68  8.9e-15
256  15  8  76  4.8e-14
512  16  10  93  6.5e-14
1024  18  13  121  1.1e-13
2048  22  19  187  3.2e-13
4096  23  22  224  1.3e-13

Table 3.6

For well conditioned matrices the number of iterations as well as the maximum displacement rank is independent of n. For the matrix 4, the number of iterations as well as the sum of the displacement ranks seems to grow as log_{2} n.

Acknowledgments. This work has been supported by the MURST project “Advanced numerical methods for scientific computing” and by the Coordinated CNR project “Toeplitz-like matrices: structures algorithms and applications”.

REFERENCES

[1] P. Amodio and F. Mazzia. Backward error analysis of cyclic reduction for the solution of tridiagonal systems. Math. Comp., 62:601–617, 1994.

[2] D. A. Bini. Using FFT-based techniques in polynomial and matrix computations: recent advances and applications. In Proceedings of FAA98, Kuwait, 1998. To appear.

[3] D. A. Bini, S. Chakravarthy, and B. Meini. A new algorithm for the design of finite capacity service units. In B. Plateau, W. J. Stewart, and M. Silva, editors, Numerical Solution of Markov Chains, 3rd International Workshop, Saragoza (Spain), Sept. 1999, pages 247– 260. Prensas Unversitarias de Zaragoza, 1999.

[4] D. A. Bini, L. Gemignani, and B. Meini. Factorization of Analytic Functions by means of Koenig’s Theorem and Toeplitz Computations. Submitted for publication, 1998.

[5] D. A. 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.

[6] D. A. 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.

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

[8] D. A. Bini and B. Meini. Using displacement structure for solving non-skip-free M/G/1 type Markov chains. In A. Alfa and S. Chakravarthy, editors, Advances in Matrix Analytic Methods for Stochastic Models - Proceedings of the 2nd international conference on matrix analytic methods, pages 17–37. Notable Publications Inc, NJ, 1998.

[9] D. A. Bini and B. Meini. Approximate displacement rank and applications. In V. Olshevsky, editor, Proceedings AMS Conference ”Structured Matrices in Operator Theory, Numerical Analysis, Control, Signal and Image Processing”, Boulder. AMS, 1999. to apppear.

[10] D. A. Bini and B. Meini. Effective methods for solving banded Toeplitz systems. SIAM J. Matrix Anal. Appl., 20:700–719, 1999.

[11] D. A. Bini and B. Meini. Solving block banded block Toeplitz systems with banded Toeplitz blocks. In F. T. Luk, editor, Advanced Signal Processing Algorithms, Architectures, and Implementations IX, volume 3807, pages 300–311. Proceedings of SPIE, 1999.

[12] D. A. Bini, B. Meini, and V. Ramaswami. Analyzing M/G/1 paradigms through QBDs: the role of the block structure in computing the matrix G. submitted for publication, 1999.

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

[14] R. R. Bitmead and B. D. O. Anderson. Asymptotically fast solution of Toeplitz and related systems of linear equations. Linear Algebra Appl., 34:103–116, 1980.

[15] H. R. Gail, S. L. Hantler, and B. A. Taylor. Solutions of the basic matrix equation for M/G/1 and G/M/1 type Markov chains. Commun. Statist. Stochastic Models, 10:1–43, 1994.

[16] 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.

[17] G. Heinig and K. Rost. Algebraic Methods for Toeplitz-like Matrices and Operators. Akademie-Verlag, Berlin, and Birkhäuser, Boston, 1984.

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

[19] A. S. Householder. The Numerical Treatment of a Single Nonlinear Equation. International Series in Pure and Applied Mathematics. McGraw-Hill Inc., New York, NY, 1970.

[20] T. Kailath. Displacement structure and array algorithms. In T. Kailath and A. H. Sayed, editors, Fast reliable algorithms for matrices with structure, chapter 1, pages 1–56. SIAM, Philadelphia, 1999.

[21] T. Kailath, S. Kung, and M. Morf. Displacement rank of a matrix. Bulletin of the American Mathematical Society, 1:769–773, 1979.

[22] T. Kailath, S. Kung, and M. Morf. Displacement ranks of matrices and linear equations. J. Math. Anal. Appl., 68:395–407, 1979.

[23] P. Lancaster and M. Tismenetski. Matrix Polynomials. Academic Press Inc., New York, 1982.

[24] G. Latouche and V. Ramaswami. A logarithmic reduction algorithm for Quasi-Birth-Death processes. J. Appl. Probability, 30:650–674, 1993.

[25] B. Meini. Efficient computation of the extreme solutions of X +A^{∗}X^{−1}A = Q, X −A^{∗}X^{−1}A = Q. Submitted for publication, 1999.

[26] M. F. Neuts. Matrix-Geometric Solutions in Stochastic Models: An Algorithmic Approach. The Johns Hopkins University Press, Baltimore, MD, 1981.

[27] V. Y. Pan, S. Branham, R. E. Rosholt, and A.-L. Zheng. Newton’s iteration for structured matrices. In T. Kailath and A. H. Sayed, editors, Fast reliable algorithms for matrices with structure, chapter 7, pages 189–210. SIAM, Philadelphia, 1999.

[28] V. Ramaswami. The generality of Quasi Birth-and-Death processes. In A. Alfa and S. Chakravarthy, editors, Advances in Matrix Analytic Methods for Stochastic Models - Proceedings of the 2nd international conference on matrix analytic methods, pages 93– 113. Notable Publications Inc, NJ, 1998.

[29] P. Tilli. Locally Toeplitz sequences: spectral properties and applications. Linear Algebra Appl., 278:91–120, 1998.

[30] P. Tilli. Asymptotic spectral distribution of Toeplitz-related matrices. In T. Kailath and A. H. Sayed, editors, Fast reliable algorithms for matrices with structure, chapter 6, pages 153–187. SIAM, Philadelphia, 1999.

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