Inverting Block Toeplitz Matrices in Block Hessenberg Form by Means of Displacement Operators: Application to Queueing Problems

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

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

Abstract

The concept of displacement rank is used to devise an algorithm for the inversion of an n × n block Toeplitz matrix in block Hessenberg form H_{n} having m × m block entries. This kind of matrices arises in many important problems in queueing theory. We explicitly relate the first and the last block rows and block columns of H_{n}^{−1} with the corresponding ones of H−1_{\frac{n}{2}}. These block vectors fully define all the entries of H_{n}^{−1} by means of a Gohberg-Semencul-like formula. In this way we obtain a doubling algorithm for the computation of H_{2^{i}}^{−1}, i = 0, 1, …, q, n = 2^{q}, where at each stage of the doubling procedure only few convolutions of block vectors must be computed. The overall cost of this computation is O(m^{2}n log n + m^{3}n) arithmetic operations with a moderate overhead constant. The same technique can be used for solving the linear system H_{n}x = b within the same computational cost. The case where H_{n} is in addition a scalar Toeplitz matrix is analyzed as well. An application to queueing problems is presented and comparisons with existing algorithms are performed showing the higher efficiency and reliability of this approach.

Key words: Block Toeplitz matrices, block Hessenberg, displacement rank, queueing problems, divide-and-conquer.

1 Introduction

Block Toeplitz matrices are frequently encountered in the modeling of problems of the real world where some sort of shift invariance holds. The shift invariance property is captured by the block Toeplitz structure where the block entries T_{i,j} of the matrix T are invariant under shift along the diagonal direction, that is T_{i,j} is a function of j − i.

In queueing theory block Toeplitz matrices play a relevant role, in fact they model a wide and important class of problems (the one described by M/G/1 type or by G/M/1 type Markov chains [14], [13]) where the block Toeplitz matrix has also the block Hessenberg structure.

One of the major problems for M/G/1 Markov chains consists in computing the first block row of the inverse of the n × n block matrix

H_{n} = [ \begin{matrix}A_{1} & A_{0} & & & ◯ \\ A_{2} & A_{1} & A_{0} \\ ⋮ & ... & ... & ... \\ A_{n−1} & ... & ... & A_{1} & A_{0} \\ A_{n} & A_{n−1} & … & A_{2} & A_{1} \end{matrix} ]  (1)

where A_{i} are m × m blocks and n generally takes very large values. In fact, for Markov chains endowed with an infinite set of states, the block Toeplitz matrix involved is infinite and we are interested in the asymptotic behavior of the inverse of H_{n}, for n → +∞.

Motivated by the large size of the matrices arising in the applications, we are interested in devising fast and efficient algorithms for inverting block Toeplitz matrices in block Hessenberg form.

In [15], [12] G. W. Stewart and G. Latouche have introduced, analyzed and implemented a recursive algorithm, that relies on the Sherman Woodbury Morrison formula, for solving systems in block Hessenberg form. This method, based on a doubling technique, has been adapted to block Toeplitz matrices and used for solving M/G/1 type Markov chains. The resulting cost of the Latouche and Stewart algorithm for solving the system H_{n}x = b is reduced to O(m^{3}n log^{2} n) arithmetic operations if Fast Fourier Transform (FFT) is used for performing computations among block Toeplitz matrices. However, this adaptation does not fully exploit the block Toeplitz structure of the matrix H_{n}. In fact, this structure is only used for the fast computation of the products between Toeplitz matrices and vectors.

The concept of displacement rank, introduced by Kailath et al. in [11], is fundamental in devising and analyzing algorithms related to Toeplitz matrices. An excellent survey on the properties and the applications of the displacement rank, together with a very rich bibliography, is presented by Kailath and Sayed in [10]. Applications to parallel computations, polynomial computations and Toeplitz/Hankel-like matrices can be found in [6], [5].

The invariance under inversion of the displacement rank of a block Toeplitz matrix allows us to fully exploit the structure of H_{n}.

In fact, by using the block displacement operator Δ(H) = ZH − HZ, where Z is the block down-shift matrix, we provide an inversion formula for the matrix H_{n} based on the properties of Δ(H_{n}). More specifically, we show that the inverse of H_{n} can be explicitly represented in terms of its first and last block rows R_{n}^{(1)}, R_{n}^{(2)}, and its first and last block columns C_{n}^{(1)}, C_{n}^{(2)}, by means of the relation

H_{n}^{−1} = L_{n}(C_{n}^{(1)})L_{n}^{T} (E_{n}^{(1)} + Z_{n}R_{n}^{(1)T} A_{0}^{T}) − L_{n}(C_{n}^{(2)})L_{n}^{T} (Z_{n}R_{n}^{(2)T} A_{0}^{T}) (2)

where L_{n}(W) denotes the n × n block lower triangular block Toeplitz matrix defined by its first block column W, and E_{n}^{(1)} denotes the nk × k matrix made up by the first k columns of the identity matrix of dimension nk.

Then we explicitly relate the block vectors R_{n}^{(1)}, R_{n}^{(2)}, C_{n}^{(1)}, C_{n}^{(2)}, to the block vectors R(1)_{\frac{n}{2}}, R(2)_{\frac{n}{2}}, C(1)_{\frac{n}{2}}, C(2)_{\frac{n}{2}}, defining H−1_{\frac{n}{2}}, and exploit these relations in order to devise a doubling algorithm for the computation of H_{n}^{−1}.

Since the equations relating block rows and block columns of H_{n}^{−1} and H−1_{\frac{n}{2}} involve products between block Toeplitz matrices and block vectors, we may implement these relations by using FFT-based techniques. In this way we devise an algorithm for the computation of H_{n}^{−1} that requires O(m^{2}n log n + m^{3}n) arithmetic operations, where the hidden constant in the asymptotic estimate has a moderate size.

Once the representation (2) of H_{n}^{−1} has been computed, we may solve any linear system H_{n}x = b with an additional cost of O(m^{2}n log n) arithmetic operations.

Moreover, we consider the important case where the matrix H_{n} is also Toeplitz in the scalar sense. This situation occurs in the modeling of certain queueing problems [9]. We show that, under this further assumption, the doubling algorithm can be strongly simplified since only the first and last rows and columns of H_{n}^{−1} must be computed. In this way we reduce the computational cost to O(m^{2} log n + mn log(mn)) arithmetic operations.

We have implemented our algorithm in Fortran 90 and compared it to the algorithm of [15] on a problem arising from the modeling of metropolitan networks [1]. Our algorithm, besides being faster than the one of [15] by a factor of 6.8, requires less storage, and therefore allows us to deal with larger values of n, and provides more accurate results.

The paper is organized as follows. In section 2 we recall the equations, based on the Sherman Woodbury Morrison formula, relating the matrices H_{n}^{−1} and H−1_{\frac{n}{2}}, on which the algorithm of [15], [12] relies. In section 3 we introduce the block displacement operator Δ(H), analyze its properties and obtain the representation formula (2) expressing the matrix H_{n}^{−1} by means of the block vectors R_{n}^{(1)}, R_{n}^{(2)}, C_{n}^{(1)}, C_{n}^{(2)}. In section 4, after recalling FFT-based procedures for performing computations between block Toeplitz matrices, we describe and analyze our algorithm. In section 5 we show an application to the solution of queueing problems, analyze the modifications needed to deal with the case of scalar Toeplitz matrices and present numerical results and comparisons.

2 A divide and conquer algorithm

Let A_{i}, i = 0, 1, …, n, be m×m matrices and consider the block Toeplitz matrix in block Hessenberg form H_{n} of (1). Let us assume for simplicity that n = 2^{q}, for a positive integer q. Suppose that det H_{2^{j}} ≠ 0, j = 0, 1, …, q, and partition the matrix H_{n} in the following way

H_{n} = [ \begin{matrix}H\frac{n}{2} & U\frac{n}{2} V T_{\frac{n}{2}} \\ T\frac{n}{2} & H\frac{n}{2} \end{matrix} ]  (3)

where

U\frac{n}{2} = [ O ; ⋮ ; O ; I ] , V\frac{n}{2} = [ A_{0}^{T} ; O ; ⋮ ; O ] , T\frac{n}{2} = [ \begin{matrix}A_{\frac{n}{2} +1} & A\frac{n}{2} & … & A_{2} \\ A_{\frac{n}{2} +2} & A_{\frac{n}{2} +1} & ... & A_{3} \\ ⋮ & ... & ... & ⋮ \\ A_{n} & A_{n−1} & … & A_{\frac{n}{2} +1} \end{matrix} ] (4)

and I denotes the m × m identity matrix. By applying the Sherman Woodbury Morrison formula [8] to the decomposition

\begin{matrix}H_{n} = [ \begin{matrix}H\frac{n}{2} & O \\ T\frac{n}{2} & H\frac{n}{2} \end{matrix} ] + [ \begin{matrix}O & U\frac{n}{2} V T_{\frac{n}{2}} \\ O & O \end{matrix} ] = S_{n} + M_{n}N_{n}^{T} , \\ M_{n} = [ \begin{matrix}U\frac{n}{2} \\ O \end{matrix} ] , N_{n} = [ \begin{matrix}O \\ V\frac{n}{2} \end{matrix} ] \end{matrix}  (5)

we immediately find the following expression for the matrix inverse of H_{n}:

H_{n}^{−1} = S_{n}^{−1} − S_{n}^{−1}M_{n}(I + N_{n}^{T} S_{n}^{−1}M_{n})^{−1}N_{n}^{T} S_{n}^{−1},  (6)

where

S_{n}^{−1} = [ \begin{matrix}H−1_{\frac{n}{2}} & O \\ −H−1_{\frac{n}{2}}T\frac{n}{2} H−1_{\frac{n}{2}} & H−1_{\frac{n}{2}} \end{matrix} ] .  (7)

The above formulae have been used by Stewart [15] in order to devise an efficient doubling method for solving the block Hessenberg system H_{n}x = b, where x and b are mn-dimensional vectors, and applied by Latouche and Stewart [12] to solving queueing problems.

In fact, partitioning x and b in two nm/2 dimensional vectors such that x^{T} = [x_{1}^{T} , x_{2}^{T}], b^{T} = [b_{1}^{T} , b_{2}^{T}], and denoting y_{1} = H−1_{\frac{n}{2}}b_{1}, y_{2} = H−1_{\frac{n}{2}}(b_{2} − T\frac{n}{2} y_{1}), Y_{3} = H−1_{\frac{n}{2}}U\frac{n}{2} , Y_{4} = −H−1_{\frac{n}{2}}T\frac{n}{2} Y_{3}, from (6), (7) we easily find that

x = [ \begin{matrix}y & 1 \\ y & 2 \end{matrix} ] − [ \begin{matrix}Y_{3} \\ Y_{4} \end{matrix} ] (I + V T_{\frac{n}{2}} Y_{4})^{−1}V T_{\frac{n}{2}} y_{2}.  (8)

In this way the solution of an n × n block Hessenberg block Toeplitz system is reduced to solving the two \frac{n}{2} × \frac{n}{2} block Hessenberg block Toeplitz systems yielding y_{1} and y_{2}, to solving the block Hessenberg block Toeplitz systems H\frac{n}{2} Y_{3} = U\frac{n}{2} and H\frac{n}{2} Y_{4} = −T\frac{n}{2} Y_{3}. Observe that all the above systems are associated with the same matrix H\frac{n}{2} and that the latter two block systems do not depend on the vector b.

Now, denote with p_{n} the cost of computing Y_{3} and Y_{4}, with s_{n} the arithmetic cost of computing x by means of (8) once the block vectors Y_{3} and Y_{4} have been computed, finally with c_{n} the overall cost of computing x. In this way from the equations (6), (8) and (7) we have

c_{n} = s_{n} + p_{n} + O(m^{2}n + m^{3})
s_{n} = 2s\frac{n}{2} + O(m^{2}n log n)
p_{n} = 2ms\frac{n}{2} + p\frac{n}{2} + O(m^{2}n log n + m^{3}n),

where the products between block Toeplitz matrices and (block) vectors are computed by means of FFT’s, as shown in section 4, for the cost O(m^{2}n log n + m^{3}n).

From the above relations we obtain s_{n} = O(m^{2}n log^{2} n), whence p_{n} = p\frac{n}{2} + O(m^{3}n log^{2} n). Therefore we deduce that p_{n} = O(m^{3}n log^{2} n) and c_{n} = O(m^{3}n log^{2} n).

In the following section, by using a suitable displacement operator, we show how the above doubling technique can be modified in order to obtain a simpler algorithm having a computational cost of O(m^{3}n + m^{2}n log n) operations.

3 Displacement properties

Let us recall some properties of displacement operators which are particularly suitable for dealing with the inversion of the matrix H_{n}. For more details on such operators we refer the reader to the pioneering paper [11], to the excellent survey given by Kailath and Sayed [10], and more specifically, concerning the operator that we use and adjust to this application, to the book [6] and to [2], [5].

Define the n × n block down-shift matrix

Z_{n} = [ \begin{matrix}O \\ I & O \\ & ... & ... \\ & & I & O \end{matrix} ] ,

where the blocks have dimension m, and consider the following block displacement operator

Δ(H) = Z_{n}H − HZ_{n}

defined for any n × n block matrix H.

Let us also introduce the following notations: the block vectors E_{n}^{(1)} = [I, O, …, O]^{T}, E_{n}^{(2)} = [O, …, O, I]^{T} denote mn × m matrices made up by the first and the last m columns of the identity matrix of order mn. Throughout we call n-dimensional block column vector any mn × m matrix, similarly we call n-dimensional block row vector any m × mn matrix. Moreover we say that the block matrix H has block displacement rank r if r is the minimum integer such that there exist block column vectors U^{(i)}, and block row vectors V^{(i)}, i = 1, …, r, satisfying the equation Δ(H) = ∑_{i=1}^{r} U^{(i)}V^{(i)}.

For m = 1 the block displacement rank of H is reduced to the well-known concept of displacement rank that coincides with the rank of Δ(H), where Δ is the (scalar) displacement operator.

It is immediate to observe that for the matrix H_{n} of (1) it holds

Δ(H_{n}) = [ I ; O ; ⋮ ; O ] [A_{0}, O, …, O] − [ O ; ⋮ ; O ; I ] [O, …, O, A_{0}].  (9)

Observe that Δ(H_{n}) only depends on the block entry A_{0} of the matrix H_{n} and is strongly structured. This fact in principle might seem a draw-back of this operator since the informations contained in the remaining blocks A_{i} are lost. However, it can be used, in the light of theorem 1, to give a simple representation of H_{n}^{−1}, well suited for block Hessenberg matrices.

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

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

The following result can be easily proved by extending to block matrices the same proof given in [6].

Theorem 1 Let K_{n} be an n×n block matrix such that Δ(K_{n}) = ∑_{i=1}^{r} U^{(i)}V^{(i)} where U^{(i)} and V^{(i)} are n-dimensional block column and block row vectors, respectively. Then we have

K_{n} = L_{n}(K_{n}E_{n}^{(1)}) − ∑_{i=1}^{r} L_{n}(U^{(i)})L_{n}^{T} (Z_{n}V^{(i)T})

where L_{n}(W) denotes the n × n block lower triangular block Toeplitz matrix defined by its first block column W.

The above result allows one to represent any matrix K_{n} as a sum of products of block lower and upper triangular block Toeplitz matrices defined by the first block column of K_{n} and by the block vectors U^{(i)}, V^{(i)} associated with the block displacement of K_{n}.

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

K_{n}^{−1} = L_{n}(K_{n}^{−1}E_{n}^{(1)}) + ∑_{i=1}^{r} L_{n}(K_{n}^{−1}U^{(i)})L_{n}^{T} (Z_{n}K_{n}^{−T} V^{(i)T}).  (11)

For a block Toeplitz matrix in block Hessenberg form the equation (11) can be much simplified. In fact, due to the special structure (9) of Δ(H_{n}) we may easily prove that

H_{n}^{−1} = L_{n}(C_{n}^{(1)})L_{n}^{T} (E_{n}^{(1)} + Z_{n}R_{n}^{(1)T} A_{0}^{T}) − L_{n}(C_{n}^{(2)})L_{n}^{T} (Z_{n}R_{n}^{(2)T} A_{0}^{T}) (12)

where C_{n}^{(1)}, C_{n}^{(2)} denote the first and the last block column, respectively, of H_{n}^{−1}, while R_{n}^{(1)}, R_{n}^{(2)} denote the first and the last block row, respectively, of the matrix H_{n}^{−1}.

In this way the inverse of H_{n} is explicitly determined as a sum of products of block triangular block Toeplitz matrices defined by the block A_{0}, by the block rows R_{n}^{(1)}, R_{n}^{(2)}, and the block columns C_{n}^{(1)}, C_{n}^{(2)}.

4 Inverting the matrix H_{n}

The results reported in sections 2 and 3 can be merged together in order to derive a doubling algorithm for the inversion of H_{n}. More precisely we determine explicit relations between the block vectors C_{n}^{(1)}, C_{n}^{(2)}, R_{n}^{(1)}, R_{n}^{(2)} defining the matrix H_{n}^{−1} and the block vectors C(1)_{\frac{n}{2}}, C(2)_{\frac{n}{2}}, R(1)_{\frac{n}{2}}, R(2)_{\frac{n}{2}} defining the matrix H−1_{\frac{n}{2}}. Moreover, since such relations only involve operations between block Toeplitz matrices and block vectors, we may devise an efficient scheme for their implementation based on performing few FFT’s. In this way we arrive at an algorithm for the computation of the inverses of H_{2^{i}}, i = 0, 1, …, q, requiring O(m^{3}n + m^{2}n log n) arithmetic operations, for n = 2^{q}, and a low storage cost; in fact, only four auxiliary block vectors need to be allocated in order to carry out our algorithm. Another interesting feature is that in the PRAM model of parallel computation the algorithm can be carried out in O(log^{2} n) parallel steps using O(m^{3}n) processors.

The following result holds

Theorem 2 The first and the last block columns C_{n}^{(1)}, C_{n}^{(2)}, and the first and the last block rows R_{n}^{(1)}, R_{n}^{(2)} of H_{n}^{−1} satisfy the following relations

\begin{matrix}C_{n}^{(1)} = [ \begin{matrix}C(1)_{\frac{n}{2}} \\ −H−1_{\frac{n}{2}}T\frac{n}{2} C(1)_{\frac{n}{2}} \end{matrix} ] − [ \begin{matrix}C(2)_{\frac{n}{2}} \\ −H−1_{\frac{n}{2}}T\frac{n}{2} C(2)_{\frac{n}{2}} \end{matrix} ] P R(1)_{\frac{n}{2}}T\frac{n}{2} C(1)_{\frac{n}{2}} \\ C_{n}^{(2)} = [ \begin{matrix}O \\ C(2)_{\frac{n}{2}} \end{matrix} ] + [ \begin{matrix}C(2)_{\frac{n}{2}} \\ −H−1_{\frac{n}{2}}T\frac{n}{2} C(2)_{\frac{n}{2}} \end{matrix} ] P R(1)_{\frac{n}{2}}E(2)_{\frac{n}{2}} \\ R_{n}^{(1)} = A_{0}[R(1)_{\frac{n}{2}}, O] − A_{0}R(1)_{\frac{n}{2}}E(2)_{\frac{n}{2}}P [R(1)_{\frac{n}{2}}T\frac{n}{2} H−1_{\frac{n}{2}}, −R(1)_{\frac{n}{2}}] \\ R_{n}^{(2)} = A_{0}[−R(2)_{\frac{n}{2}}T\frac{n}{2} H−1_{\frac{n}{2}}, R(2)_{\frac{n}{2}}] + A_{0}R(2)_{\frac{n}{2}}T\frac{n}{2} C(2)_{\frac{n}{2}}P [R(1)_{\frac{n}{2}}T\frac{n}{2} H−1_{\frac{n}{2}}, −R(1)_{\frac{n}{2}}], \end{matrix} (13)

where P = (I + A_{0}R(1)_{\frac{n}{2}}T\frac{n}{2} C(2)_{\frac{n}{2}})^{−1}A_{0}.

Proof. The formulae follow directly from (5–7). □

In order to describe the algorithm for the computation of H_{n}^{−1}, represented in terms of the block vectors R_{n}^{(1)}, R_{n}^{(2)}, C_{n}^{(1)}, C_{n}^{(2)}, based on Theorem 2, we need to recall some procedures for performing computations involving block Toeplitz matrices. For this purpose we denote with y = (y_{0}, …, y_{n−1}) = DFT(x) the Discrete Fourier Transform of order n of the vector x = (x_{0}, …, x_{n−1}) defined by y_{i} = ∑_{j=0}^{n−1} ω^{ij}x_{j}, i = 0, …, n − 1, where ω = cos\frac{2π}{n} + i sin\frac{2π}{n}, and i is the imaginary unit such that i^{2} = −1. Analogously x = IDFT(y) denotes the Inverse Discrete Fourier Transform of order n such that x_{i} = \frac{1}{n} ∑_{j=0}^{n−1} ω^{−ij}y_{j}, i = 0, …, n − 1. We recall that the computation of a DFT of order n = 2^{q} of a real vector can be performed by means of the base 2 FFT algorithm in about \frac{5}{2}n log_{2} n arithmetic operations between real numbers [7]. Similarly, the cost of an IDFT is about \frac{5}{2}n log_{2} n + n real arithmetic operations in the case where the transformed vector is real.

The above notations can be extended to block vectors. Given the block vector X = (X_{0}, …, X_{n−1}), where X_{s} = (x_{i,j}^{(s)}) are m × m matrices, we define Y = (Y_{0}, …, Y_{n−1}) = DFT(X) the block vector such that Y_{s} = (y_{i,j}^{(s)}) are m×m matrices and (y_{i,j}^{(0)}, …, y_{i,j}^{(n−1)}) = DFT(x_{i,j}^{(0)}, …, x_{i,j}^{(n−1)}), i, j = 1, …, m. Analogously we set X = IDFT(Y) if (x_{i,j}^{(0)}, …, x_{i,j}^{(n−1)}) = IDFT(y_{i,j}^{(0)}, …, y_{i,j}^{(n−1)}). In this way the computation of DFT(X) and IDFT(Y) is reduced to the computation of m^{2} DFT’s and IDFT’s, respectively, for the cost of about \frac{5}{2}m^{2}n log_{2} n and \frac{5}{2}m^{2}n log_{2} n + mn real operations, respectively.

Remark 3 Observe that the product of a block Toeplitz matrix and a vector can be reduced to computing the product of two matrix polynomials, that is, polynomials whose coefficients are matrices. In fact, if P (z) = ∑_{i=0}^{n−1} P_{i}z^{i}, Q(z) = ∑_{i=0}^{n/2−1} Q_{i}z^{i}, are matrix polynomials, where P_{i}, Q_{i} are m×m matrices, then, by comparing the terms of degree n/2, …, n − 1 in the equation R(z) = P (z)Q(z) mod (z^{n} − 1), we arrive at the following equation involving a block Toeplitz matrix

[ R_{n/2} ; ⋮ ; R_{n−1} ] = [ \begin{matrix}P_{n/2} & … & P_{1} \\ ⋮ & ... & ⋮ \\ P_{n−1} & … & P_{n/2} \end{matrix} ] [ Q_{0} ; ⋮ ; Q_{n/2−1} ] .

The product of two matrix polynomials can be computed by the following algorithm.

Algorithm 1. Computation of matrix polynomial product modulo z^{n} − 1.

Input. Positive integers q, m and the m × m matrices P_{0}, …, P_{n−1}, Q_{0}, …, Q_{n−1}, coefficients of the matrix polynomials P (z) = ∑_{i=0}^{n−1} P_{i}z^{i},

Q(z) = ∑_{i=0}^{n−1} Q_{i}z^{i}, where n = 2^{q}.

Output. The m × m matrices R_{0}, …, R_{n−1}, coefficients of the matrix polynomial R(z) = ∑_{i=0}^{n−1} R_{i}z^{i} such that R(z) = P (z)Q(z) mod (z^{n} − 1).

Computation.

1. (Evaluation) Compute the entries of the 2n matrices U_{s} = P (ω^{s}), V_{s} = Q(ω^{s}), s = 0, …, n − 1, ω = cos\frac{2π}{n} + i sin\frac{2π}{n}, in the following way

\begin{matrix}(U_{0}, …, U_{n−1}) = DFT(P_{0}, …, P_{n−1}), \\ (V_{0}, …, V_{n−1}) = DFT(Q_{0}, …, Q_{n−1}). \end{matrix}

2. Compute the n matrix products W_{s} = U_{s}V_{s}, s = 0, …, n − 1, such that W_{s} = P (ω^{s})Q(ω^{s}) = R(ω^{s}).

3. (Interpolation) Compute the entries of the matrices R_{0}, …, R_{n−1} by means of the equation (R_{0}, …, R_{n−1}) = IDFT(W_{0}, …, W_{n−1}).

The cost of Algorithm 1 is O(m^{3}n + m^{2}n log n) ops. More precisely, for real input matrices the cost is less than \frac{5}{2}m^{2}n log_{2} n + n(3m^{3} + 2m^{2}) (compare [4], [3]).

Remark 4 By using equation (12) and Algorithm 1 it is possible to compute the product of the matrix H_{n}^{−1} and a block column vector by means of 10m^{2} FFT’s of order 2n and about 4n complex m × m matrix products. This algorithm can be easily adjusted to the computation of the product of a block row vector and H_{n}^{−1}. Moreover, for the computation of the products of H_{n}^{−1} with k block row and/or column vectors the number of FFT’s can be reduced to 6k + 4 while the number of complex matrix products is 4kn. In fact, the DFT’s of the four block vectors defining H_{n}^{−1} can be computed once for all.

Now we are ready to describe the main algorithm of this paper.

Algorithm 2. Computation of the inverse of a block Hessenberg block Toeplitz matrix H_{n}.

Input. A positive integer q and the m × m blocks A_{i}, i = 0, 1, …, n, n = 2^{q}, such that det H_{2^{j}} ≠ 0, j = 0, 1, …, q.

Output. The block vectors C_{2^{j}}^{(1)}, C_{2^{j}}^{(2)}, R_{2^{j}}^{(1)}, R_{2^{j}}^{(2)}, defining the first and the last block columns and block rows, respectively, of the matrix H_{2^{j}}^{−1}, for j = 0, 1, …, q.

Computation.

1. Compute H_{1}^{−1} = A_{1}^{−1} and set R_{1}^{(1)} = R_{1}^{(2)} = C_{1}^{(1)} = C_{1}^{(2)} = A_{1}^{−1}.

2. For j = 0, 1, …, q − 1 compute C_{2^{j+1}}^{(1)}, C_{2^{j+1}}^{(2)}, R_{2^{j+1}}^{(1)}, R_{2^{j+1}}^{(2)} by means of Theorem 2, by calculating:

(a) X^{(1)} = T_{2^{j}}C_{2^{j}}^{(1)}, X^{(2)} = T_{2^{j}}C_{2^{j}}^{(2)}, Y^{(1)} = R_{2^{j}}^{(1)}T_{2^{j}}, Y^{(2)} = R_{2^{j}}^{(2)}T_{2^{j}};

(b) H_{2^{j}}^{−1}X^{(1)}, H_{2^{j}}^{−1}X^{(2)}, Y^{(1)}H_{2^{j}}^{−1}, Y^{(2)}H_{2^{j}}^{−1};

(c) P = (I + A_{0}R_{2^{j}}^{(1)}X^{(2)})^{−1}A_{0};

(d) C_{2^{j+1}}^{(1)}, C_{2^{j+1}}^{(2)}, R_{2^{j+1}}^{(1)}, R_{2^{j+1}}^{(2)} by means of (13).

For the computation of the products T_{2^{j}}C_{2^{j}}^{(1)}, T_{2^{j}}C_{2^{j}}^{(2)}, R_{2^{j}}^{(1)}T_{2^{j}}, R_{2^{j}}^{(2)}T_{2^{j}}, at stage 2a of Algorithm 2, involving the block Toeplitz matrix T_{2^{j}}, we use Algorithm 1 where the FFT of the matrix T_{2^{j}}, that appears in each product, is computed once. In this way we have to compute 9m^{2} FFT’s of order 2^{j+1} and about 2^{j+2} products of m × m complex matrices.

For the computation of the four products involving the matrix H_{2^{j}}^{−1} at stage 2b of Algorithm 2, we use again Algorithm 1 and (12) in the light of Remark 4, where the four FFT’s associated with the matrix H_{2^{j}}^{−1} are computed only once. In this way the computation is reduced to performing 12m^{2} FFT’s of order 2^{j+1} and about 10 · 2^{j} complex m × m matrix products.

Since the cost of an FFT of a real vector of order n is about \frac{5}{2}n log_{2} n real arithmetic operations and the product of two complex matrices can be performed by means of three products of real matrices, we find that the overall cost of Algorithm 2 is roughly given by 42m^{3}n + 100m^{2}n log_{2} n arithmetic operations.

Once the block vectors R_{n}^{(1)}, R_{n}^{(2)}, C_{n}^{(1)}, C_{n}^{(2)} have been computed, the solution x of the system H_{n}x = b, can be obtained by means of (12) and remark 3, by performing O(m^{2}) FFT’s of order n, thus, leaving unchanged the asymptotic cost O(m^{3}n + m^{2}n log n).

5 Application to queueing problems

Algorithm 2 proposed in Section 4 can be applied to solve nonlinear matrix equations arising in queueing problems. Specifically, it is fundamental for solving Markov chains of M/G/1 type (for details on M/G/1 type Markov chains we refer the reader to [14]) the computation of the minimal nonnegative solution G of the nonlinear matrix equation

X = ∑_{i=0}^{∞} X^{i}B_{i},  (14)

where B_{i}, i ≥ 0, are nonnegative m × m matrices such that ∑_{i=0}^{∞} B_{i} is column stochastic, and X is an m × m matrix. In [12] it is observed that the seeked solution G satisfies the following equation involving an infinite block Hessenberg block Toeplitz matrix

[G, G^{2}, G^{3}, …] [ \begin{matrix}I − B_{1} & −B_{0} & & & ◯ \\ −B_{2} & I − B_{1} & −B_{0} \\ −B_{3} & −B_{2} & I − B_{1} & −B_{0} \\ ⋮ & ... & ... & ... & ... \end{matrix} ] = [B_{0}, O, O, …].

(15)

Moreover, if we consider the solution of the block n × n truncated system

\begin{matrix}[X_{1}^{(n)}, X_{2}^{(n)}, …, X_{n}^{(n)}]H_{n} = [B_{0}, O, …, O] \\ H_{n} = [ \begin{matrix}I − B_{1} & −B_{0} & & & ◯ \\ −B_{2} & I − B_{1} & −B_{0} \\ ⋮ & ... & ... & ... \\ −B_{n−1} & ... & ... & I − B_{1} & −B_{0} \\ −B_{n} & −B_{n−1} & … & −B_{2} & I − B_{1} \end{matrix} ] \end{matrix}  (16)

then {X_{1}^{(n)}}_{n} is a sequence of nonnegative matrices that monotonically converges to G [12]. In order to solve the infinite system (15) it is suggested to solve the above systems, for n equal to 2, 4, …, 2^{q} for a sufficiently large q.

Since the solution of (16) is given by the first block row of H_{n}^{−1} multiplied on the left by B_{0}, we may apply Algorithm 2 of section 4 for the solution of the matrix equation (14). The cost of this computation amounts to O(m^{3}n + m^{2}n log n) real arithmetic operations.

In certain important applications [9] H_{n} is a scalar Toeplitz matrix, consequently, B_{0} is a lower triangular Toeplitz matrix that we may assume nonsingular (its non-singularity depends on the choice of the block size m). In this case it has been proved by Gail et al. [9], that the matrix G can be represented as G = F^{m}, where F is a suitable Frobenius matrix. In particular, G is such that rank(Δ(G)) = 1, that is, its displacement rank is 1, according to the operator Δ(A) = ZA − AZ, where Z is the scalar down-shift matrix.

Here we show that, for a scalar Toeplitz matrix H_{n} in block Hessenberg form, Algorithm 2 can be strongly simplified. In fact, since H_{n} is Toeplitz we may write H_{n}^{−1} as H_{n}^{−1} = L(u_{n}^{(1)})L^{T} (v_{n}^{(1)}) − L(u_{n}^{(2)})L^{T} (v_{n}^{(2)}) for suitable mn dimensional vectors u_{n}^{(i)}, v_{n}^{(i)}, i = 1, 2 (say, by using the Gohberg Semencul formula [10], or the scalar version of the formula (12)) and L(u) denotes the lower triangular Toeplitz matrix defined by its first column u..

In this way it is sufficient to relate the vectors u_{n}^{(i)}, v_{n}^{(i)}, with u(i)_{\frac{n}{2}}, v(i)_{\frac{n}{2}}, for i = 1, 2. This relation is implicitly provided by (5–7). We observe that the computation of u_{n}^{(i)}, v_{n}^{(i)} is reduced to performing a finite number of products of scalar Toeplitz matrices and vectors of dimension nm except for the computation of the inverse of I + N_{n}^{T} S_{n}^{−1}M_{n}. For this purpose we have the following result:

Theorem 5 The displacement rank of I + N_{n}^{T} S_{n}^{−1}M_{n} is at most 6, i.e., there exist vectors u^{(i)}, v^{(i)} such that Δ(I + N_{n}^{T} S_{n}^{−1}M_{n}) = ∑_{i=1}^{6} u^{(i)}v^{(i)T}.

Proof. It is sufficient to compute the rank of Δ(I + N_{n}^{T} S_{n}^{−1}M_{n}), i.e., rank(Δ(N_{n}^{T} S_{n}^{−1}M_{n})). Due to the structure of N_{n} and M_{n}, we have N_{n}^{T} S_{n}^{−1}M_{n} = A_{0}(S_{n}^{−1})_{\frac{n}{2} +1, \frac{n}{2}} , where (S_{n}^{−1})_{\frac{n}{2} +1, \frac{n}{2}} is the (\frac{n}{2} + 1, \frac{n}{2}) block entry of S_{n}^{−1}. Now, from the definition of displacement rank it follows that rankΔ((S_{n}^{−1})_{\frac{n}{2} +1, \frac{n}{2}} ) ≤ rankΔ(S_{n}^{−1})+2, moreover, since S_{n} is a 2×2 block triangular block Toeplitz matrix with Toeplitz blocks, we have rankΔ(S_{n}) ≤ 4, and therefore rankΔ(S_{n}^{−1}) ≤ 4. Whence we deduce that rankΔ((S_{n}^{−1})_{\frac{n}{2} +1, \frac{n}{2}} ) ≤ 6. Since A_{0} is a lower triangular Toeplitz matrix then it has null displacement rank and from the relation Δ(AB) = AΔ(B) + Δ(A)B, that can be easily verified, it follows rankΔ(A_{0}(S_{n}^{−1})_{\frac{n}{2} +1, \frac{n}{2}} ) ≤ rankΔ((S_{n}^{−1})_{\frac{n}{2} +1, \frac{n}{2}} ) ≤ 6. □

From the above theorem and from the equations (6–7) applied in the scalar case, it follows that the inverse of I+N_{n}^{T} S_{n}^{−1}M_{n} can be represented as the sum of at most 7 products of lower and upper triangular Toeplitz matrices respectively. Thus its inverse can be computed by means of fast algorithms which require O(m^{2}) arithmetic operations [10]. In this way, due to the above properties, each doubling step of size 2^{j} of Algorithm 2 applied to a scalar Toeplitz matrix is performed with a finite number of FFT’s of order m2^{j} and with the inversion of a matrix having displacement rank at most 4. Whence this computation can be carried out at the arithmetic cost O(m2^{j} log(m2^{j})) + O(m^{2}), and the overall cost is O(m^{2} log n) + O(mn log(mn)).

Algorithm 2 has been implemented in Fortran 90 and compared with the Stewart implementation of the algorithm given in [15], on a problem arising from the modeling of a metropolitan network [1], where the goal of the computation is solving the matrix equation (14). The size of the blocks in this problem is m = 16 and the maximum value of n that we have tested for our algorithm is n = 4096.

The experiments have been performed on an Alpha server. The following table reports, for the values of n ranging from 256 to 4096, the number of cpu seconds, required by our algorithm (algorithm DR: Displacement Rank) and the residual error; it reports also the time required by the Latouche–Stewart algorithm (denoted by LS) to reach the same residual error of our algorithm, and the ratio between the times. The symbol “*” denotes that we were not able to reach the required accuracy for lack of memory.

DR  LS
n time residual  time residual ratio
256 8.7 2.5 · 10^{−3} 38.3 2.5 · 10^{−3} 4.4
512 18.1 4.0 · 10^{−4} 103.5 4.0 · 10^{−4}  5.7
1024 39.2 1.7 · 10^{−5} 264.9 1.7 · 10^{−5}  6.8
2048 89.3 3.4 · 10^{−8}  *  *  *
4096 193.4 1.3 · 10^{−13}  *  *  *

Table 1: Comparison of the algorithms

Our algorithm, besides being faster than Stewart’s algorithm, requires less storage and therefore allows us to deal with larger values of n, and provides more accurate results.

References

[1] G. Anastasi, L. Lenzini, B. Meini, “Performance evaluation of a worst case model of the Metaring MAC Protocol with global fairness”, Performance Evaluation, to appear.

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

[3] D. Bini, B. Meini, “On cyclic reduction applied to a class of Toeplitz-like matrices arising in queueing problems”, Computations with Markov Chains, W. J. Stewart Ed., pp. 21–38, Kluwer Academic Publisher s, 1996.

[4] D. Bini, B. Meini, “On the solution of a nonlinear matrix equation arising in queueing problems”, SIAM J. Matrix Anal. Appl., 17 (1996), pp. 906– 926.

[5] D. Bini, V. Pan, “Improved parallel computations with Toeplitz-like and Hankel matrices”, Linear Algebra Appl., 188,189: 3–29 (1993).

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

[7] D.F. Elliott, K.R. Rao, Fast Transform Algorithms, Analyses, Applications, Academic Press, New York, 1982.

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

[9] H.R. Gail, S.L. Hantler, B.A. Taylor, “Non-skip-free M/G/1 and G/M/1 Type Markov Chains”, Adv. Appl. Prob., to appear.

[10] T. Kailath, A.H. Sayed, “Displacement structure: theory and applications”, SIAM Rev., 37 (1995), pp. 297–386.

[11] T. Kailath, A. Viera, M. Morf, “Inverses of Toeplitz operators, innovations and orthogonal polynomials”, SIAM Rev., 20, (1978), pp. 106–119.

[12] G. Latouche, G.W. Stewart, “Numerical methods for M/G/1 type queues”, Computations with Markov Chains, W. J. Stewart Ed., pp. 571– 581, Kluwer Academic Publishers, 1996.

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

[14] M.F. Neuts, Structured Stochastic Matrices of M/G/1 Type and Their Applications, Dekker Inc., New York, 1989.

[15] G. W. Stewart, “On the solution of block Hessenberg systems”, Numer. Lin. Algebra with Appl., 2 (1995), pp. 287–296.