Solving block banded block Toeplitz systems with structured blocks: new algorithms and open

problems

Dario Andrea Bini and Beatrice Meini^{∗}

Abstract

Some of the main results concerning Toeplitz matrix computations are recalled, in particular the recent techniques of the authors for the solution of banded Toeplitz systems are revisited. These tools are used for the design and analysis of new algorithms for solving block banded block Toeplitz systems with blocks that are banded Toeplitz matrices, and for solving certain block tridiagonal block Toeplitz systems with blocks having a tensor product structure. Applications to queueing theory, polynomial computations and image restoration are shown.

1 Matrix structures and their role

The mathematical modeling of problems of the real world often leads to solving discrete problems endowed with specific structures and with very large dimensions. For instance, it is very frequent to encounter very large linear systems whose matrix is defined by few parameters that characterize the matrix 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 and computer algebra.

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.

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

E-mail: bini@dm.unipi.it

Indeed, the strongest is the structure the most 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. This task is not trivial as it could seem at first. In fact, it happens that algorithms designed for the exploitation of the outer structures usually 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. Then we present some recent algorithms, based on the Toeplitz matrix technology, for the solution of certain applicative problems. In particular, we present problems with a two-level structure, like solving block banded block Toeplitz systems with banded Toeplitz blocks, or with blocks having a Kronecker 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.

1.1 Structures in queueing problems

In queueing problems modeled by 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 can be strongly structured, according to the model of the queueing problem. In particular, in the case of QBD (Quasi-Birth-Death) processes [12], P is a block tridiagonal block Toeplitz matrix, i.e,

Q = [ \begin{matrix}A & C & & ◯ \\ B & A & C \\ & B & A & ... \\ ◯ & & ... & ... \end{matrix} ] ,  (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 [12, 1], 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}β,  (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 [7, 3]; 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} ] ,

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. Thus, the matrix Q can be reblocked into M × M block matrices, where M = max(k, N − k). The matrix obtained in this way is a block tridiagonal, block Toeplitz matrix, with structured block entries. Hence, the problem of the solution of system Qx = b can be reduced to the problem of the solution of a system with a matrix having the structure (1), where the blocks A, B, C are themselves structured.

1.2 Structures in polynomial computations

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}.

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 [9] 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(σ^{n}), 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 [2]).

Thus the problem of approximating a factor of a polynomial is reduced to performing Toeplitz computations.

1.3 Structures in image processing

In the problems of image restoration, one of the main issues is to compute the entries of the matrix X = (x_{i,j})_{i=1,...,m,j=1,...,n}, such that

b_{i,j} = ∑_{h=−M}^{M} ∑_{k=−M}^{M} α_{h,k}x_{i+h,j+k}, i = 1, …, m, j = 1, …, n,

where we assume x_{i,j} = 0 if i and j are out of range. Here B = (b_{i,j}) represents the blurred image, X is the original image and α = (α_{i,j})_{i,j=−M,...,M} is the point-spread-function (psf) that represents the blur action; more precisely, α_{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 columnwise as vectors x and b 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. The block banded structure is due to the locality of the psf.

2 Some tools for structure analysis

In this section we recall some basic tools for the design of efficient algorithm for solving certain structured systems.

2.1 Displacement rank

The concept of displacement operators and displacement rank, introduced in [11, 10], and elaborated in [8, 5], is a powerful tool for dealing with Toeplitz matrices. Here we recall the main ideas and results.

Let us consider the operator

Δ(A) = AZ − ZA, Z = (z_{i,j}), z_{i+1,i} = 1, z_{i,j} = 0, otherwise, (3)

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 defined by its first column a. It holds

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

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 (4) is called the displacement representation of A. Equivalent representations can be given in terms of different operators.

A simple consequence of the relation

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

that can be directly derived from (3), 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}).

The above representations allow us to make elementary operations among Toeplitz-like matrices at a low cost by means of FFT. In particular, computing the product of matrices with displacement rank r 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 [6].

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 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.

Definition 2.1 Let ||·|| denote the Euclidean norm. 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).

Thus, from the above definition, 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})

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

The following properties hold.

Theorem 2.2 Let Δ(A) = UΣV^{H} be the singular value decomposition of Δ(A), where σ_{1} ≥ … ≥ σ_{n} ≥ 0 are the singular values of Δ(A). Then Drk_{ϵ}(A) = r if and only if σ_{r} > ϵ ≥ σ_{r+1}. Moreover Rdrk_{ϵ}(A) = r if and only if σ_{r} > ϵ||A|| ≥ σ_{r+1}.

The following inequalities are an immediate consequence of the definition:

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

(6)

where Cond(A) = ||A|| · ||A^{−1}||.

2.3 Cyclic reduction for block tridiagonal block Toeplitz systems

Consider the system

Qx = b,  (7)

where Q is a finite block tridiagonal block Toeplitz matrix, of block size 2^{q}, say, obtained by truncating (1). 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 (7) 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} ]  (8)

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 & & & & ◯ \\ & . & . & . \\ ◯ & & & & 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 (8) 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}  (9)

Denoting

Q^{(1)} = D_{2}^{(0)} − L^{(0)}D_{1}^{(0)^{−1}}U^{(0)}  (10)

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

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

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

It is easy to observe 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 (11) 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}F^{(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)}, F^{(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)}, F^{(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)} \\ F^{(j+1)} = F^{(j)} − C^{(j)}A^{(j)^{−1}}B^{(j)}. \end{matrix}

It is useful to express this 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, ….

From (10) we may easily deduce that

\{ \begin{matrix}ϕ^{(j+1)}(z) = zA^{(j)} − (C^{(j)} + zB^{(j)})A^{(j)^{−1}}(C^{(j)} + zB^{(j)}) \\ F^{(j+1)} = F^{(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 that

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

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 (12) is useful to prove structural properties of the block entries obtained at each step of cyclic reduction. Indeed, from (5) 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.

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

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 = (q_{i,j}) the matrix of size n having entries q_{i,j} = a_{i−j}. Q is a banded Toeplitz matrix, with bandwidth m = max{m_{1}, −m_{2}}. Assuming that n = m2^{q}, the matrix Q can be partitioned into m × m blocks yielding a block tridiagonal Toeplitz matrix with the structure (1). In order to solve Qx = b, we apply the cyclic reduction algorithm to the block tridiagonal system, as described in section 2.3. It is easy to verify that the matrix polynomial ϕ(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 formal matrix power series ψ(z) = ϕ(z)^{−1} is a z-circulant matrix. The matrix power series ψ^{(j)}(z), having matrix coefficients taken from the matrix coefficients of ψ(z), are not generally z-circulant matrices. However, since z-circulant matrices are Toeplitz and therefore the matrix coefficients of ψ(z) are Toeplitz, the coefficients of ψ^{(j)}(z) are Toeplitz matrices, for any j.

This simple fact allows us to say that ψ^{(j)}(z), and therefore ϕ^{(j)}(z), has displacement rank at most 2 for any z and for any j. In this way, efficient formulae for the computation of ϕ^{(j)}(z) can be easily devised [4].

3.2 The case of PH/PH/1 queueing problems

For the matrix of (1) with blocks (2), we may prove that [1] 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}})
F^{(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 algorithm has been implemented and the results of the numerical tests proved its reliability and efficiency. More details can be found in [1].

3.3 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 power series ψ(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.2. Indeed, under the assumptions that guarantee the convergence to zero of the blocks C^{(j)}, B^{(j)}, 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)}, F̃^{(j)}, of the blocks A^{(j)}, B^{(j)}, C^{(j)}, F^{(j)}, expressed in terms of their relative ϵ-displacement representation is based on the above property.

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

Algorithm. Input: A relative ϵ-displacement representation of Ã^{(j)}, B̃^{(j)}, C̃^{(j)}, F̃^{(j)}. Output: A relative ϵ-displacement representation of Ã^{(j+1)}, B̃^{(j+1)}, C̃^{(j+1)},

F̃^{(j+1)}.

Computation:

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

2. From the latter, and from the relative ϵ-displacement representation of B̃^{(j)} compute a relative θ-displacement representation of

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

for a suitable θ, according to formulae (6).

Let h be the rank of such representation and S_{1} and S_{2} be the mM × h matrices whose columns define a representation, i.e., such that

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

3. 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}^{H}. Compute the singular value decomposition W = UΣV^{H} 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 Σ.

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

5. 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 F̃^{(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.

Acknowledgments. This work has been supported by the MURST project “Advanced numerical methods for scientific computing”.

References

[1] D. A. Bini, S. Chakravarthy, and B. Meini. A new algorithm for the design of finite capacity service units. In Proceedings of the Third International Conference on the Numerical Solution of Markov Chains, Saragoza, Spain, Sept. 1999, in printing.

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

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

[4] D. A. Bini and B. Meini. Effective methods for solving banded Toeplitz systems. SIAM J. Matrix Anal. Appl., 1999. To appear.

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

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

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

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

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

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

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

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