Efficient computation of the extreme solutions of X + A^{∗}X^{−1}A = Q and X − A^{∗}X^{−1}A = Q

Beatrice Meini^{∗}

September 7, 2000

Abstract

We propose a new quadratically convergent algorithm, having a low computational cost per step and good numerical stability properties, which allows the simultaneous approximation of the extreme solutions of the matrix equations X + A^{∗}X^{−1}A = Q and X − A^{∗}X^{−1}A = Q. The algorithm is based on the cyclic reduction method.

Key words: Matrix equation, cyclic reduction, block Toeplitz matrix. AMS Classification: 15A24, 65F10, 65H05.

1 Introduction

In this paper we are concerned with the design and analysis of a quadratically convergent algorithm for computing the extremal solutions of the matrix equations

X + A^{∗}X^{−1}A = Q  (1)

and

X − A^{∗}X^{−1}A = Q,  (2)

where Q is an m × m Hermitian positive definite matrix, A is an m × m matrix, and A^{∗} denotes the conjugate transpose of A. An algorithm is quadratically convergent (or, equivalently, has a double exponential convergence) if it generates a sequence {X_{n}}_{n} of approximations to the solution X such that ||X − X_{n}|| ≤ γσ^{2^{n}}, for positive constants σ < 1 and γ, where || · || denotes any matrix norm; if the bound ||X − X_{n}|| ≤ γσ^{n} holds, then we say that the convergence is linear.

In [10] it is shown that, if (1) has a positive definite solution X, then there exist minimal and maximal solutions X_{−} and X_{+}, respectively, such that 0 < X_{−} ≤ X ≤ X_{+} for any positive definite solution X. Here and hereafter, if X and Y are Hermitian matrices, X ≤ Y (X < Y) means that Y − X is positive semidefinite (definite).

∗Dipartimento di Matematica, Università di Pisa, via Buonarroti 2, 56127 Pisa, Italy, meini@dm.unipi.it

Concerning equation (2), in [11] it is proved that there always exists a unique positive definite solution X_{+}, which is the maximal one, and, if A is nonsingular, there exists a unique negative definite solution X_{−}, which is the minimal one. Here maximal and minimal must be intended according to the ordering for Hermitian matrices that we have introduced above.

We will refer to X_{−} and X_{+} as to the extreme solutions of (1), or of (2). Equations (1) and (2) arise in a wide variety of research areas, which include

control theory, ladder networks, dynamic programming, stochastic filtering and statistics (see [2, 28] for a list of references concerning (1), and [11] for references about equation (2)). Recently, theoretical properties of the extreme solutions and numerical methods for their computation have been investigated by several authors (see the papers [15, 14, 11, 28, 29, 9, 10, 2]). The available numerical methods are mainly based on fixed point iterations, or on applications of Newton’s algorithm.

In the recent paper [15] the authors analyze convergence properties of fixed point iterations and of Newton’s method for computing the maximal solutions of (1) and (2). They show that the rate of convergence of those algorithms is related to the spectral radius ρ(X_{+}^{−1}A) of the matrix X_{+}^{−1}A. The convergence of fixed point iterations is linear when ρ(X_{+}^{−1}A) < 1. However, if ρ(X_{+}^{−1}A) is close to 1, the convergence may be very slow and it is more efficient to apply Newton’s method, which has a double exponential convergence. On the other hand, the drawback of the application of Newton’s method is the large computational cost per iteration; in fact, at each step the solution X of the matrix equation X + BXC = D, where X, B, C, D are m × m matrices, must be computed.

In this paper we devise a new algorithm that provides the simultaneous approximation of the extreme solutions of (1), or of (2). The algorithm has a double exponential convergence, like Newton’s method, and a low computational cost for each step, like fixed point iterations. Moreover, it shows good numerical stability properties.

The idea on which this algorithm is based arises from the closeness of equations (1) and (2) to certain problems related to Markov chains. Indeed, a similar (even though more complicated) matrix equation, is encountered in the numerical solution of Markov chains described by Quasi Birth and Death processes (see [23]), where the fundamental problem is the computation of the minimal nonnegative solution G of

G = A_{0} + A_{1}G + A_{2}G^{2},  (3)

with A_{0}, A_{1}, A_{2} nonnegative matrices such that A_{0} + A_{1} + A_{2} is stochastic (i.e. (A_{0} + A_{1} + A_{2})e = e and e is the vector having all the entries equal to one). Here, minimal means that any other solution is component-wise larger than the minimal one. The matrix G is related to the minimal nonnegative solution U of the matrix equation

I − U + A_{2}(I − U)^{−1}A_{0} = I − A_{1},  (4)

in fact, G = (I − U)^{−1}A_{0} (see [16, 18]). In the field of Markov chains, several algorithms have been developed for the numerical solution of (3), see [17, 19, 5, 6, 1, 20]. One of the most efficient methods, based on cyclic reduction [13], is proposed in [5, 6] and consists in generating a sequence of nonlinear matrix equations of the form (3), such that the matrix coefficients quadratically converge to the coefficients of a linear equation that can be explicitly solved. The procedure provides the desired solution G and the corresponding matrix U.

In this paper we show that the cyclic reduction algorithm can be efficiently applied for solving also the matrix equations (1) and (2). We observe that the matrix X_{+}^{−1}A solves a quadratic matrix equation of the form (3) and, according to the ideas developed in [5, 6], we rewrite the matrix equation in terms of an infinite block tridiagonal block Toeplitz system. By applying the cyclic reduction algorithm to this system, we generate a sequence of infinite block tridiagonal block systems, which are block Toeplitz except for the block entry in position (1, 1). Due to the nice spectral properties of X_{+}^{−1}A, we show that the sequence of matrices defining the infinite systems quadratically converges to a block diagonal matrix, and the first block diagonal entry of these matrices quadratically converges to X_{+}. More precisely, if we denote by X_{n} the block entry in position (1, 1) of the infinite matrix obtained at the n-th step of cyclic reduction, we prove that for equation (1), in the case where ρ(X_{+}^{−1}A) < 1 we have ||I −X_{n}X_{+}^{−1}|| = O((σ +ϵ)^{2·2^{n}}) for any ϵ > 0 and for any matrix norm ||·||, where σ = ρ(X_{+}^{−1}A). For equation (2), it always holds ρ(X_{+}^{−1}A) < 1, thus the double exponential convergence is always guaranteed. If A is nonsingular, the algorithm provides the simultaneous approximation of X_{−}, for both equations. Each step of the algorithm only requires the solution of two linear systems and the computation of three matrix products; moreover, the conditioning of the matrices defining the linear systems is bounded from above by a constant.

The performed experiments show that our algorithm is much faster than the existing ones in many cases.

The paper is organized as follows. In Section 2 we recall theoretical results about the existence and properties of the solutions. In Sections 3 and 4 we present the new algorithm and perform some comparisons with the known ones for the solution of (1) and (2), respectively. In Section 5 we report some numerical results.

2 Solutions of the matrix equations X+A^{∗}X^{−1}A = Q and X − A^{∗}X^{−1}A = Q

In this section we recall conditions about the existence of the extreme solutions of (1) and (2), and some spectral properties of the matrix X_{+}^{−1}A, which will be used in the subsequent sections to show the convergence of our algorithm.

2.1 The equation X + A^{∗}X^{−1}A = Q

Necessary and sufficient conditions for the existence of a positive definite solution of (1) are provided in [10]. More specifically, let us introduce the rational matrix function

ψ(λ) = λA + Q + λ^{−1}A^{∗},  (5)

defined on the unit circle C of the complex plane, which is Hermitian for any λ ∈ C. This function is said regular if there exists at least a λ ∈ C such that det ψ(λ) ≠ 0.

The following fundamental results hold [10]:

Theorem 2.1 Equation (1) has a positive definite solution X if and only if ψ(λ) is regular and ψ(λ) ≥ 0 for all λ ∈ C. Moreover, if equation (1) has a positive definite solution, then it has a maximal and minimal solution X_{+} and X_{−}, respectively.

Theorem 2.2 Suppose that A is nonsingular. Then X solves

X + A^{∗}X^{−1}A = I  (6)

if and only if Y = I − X solves

Y + AY^{−1}A^{∗} = I.  (7)

In particular, if Y_{+} is the maximal solution of (7), then X_{−} = I − Y_{+} is the minimal solution of (6).

The latter theorem can be easily generalized to the matrix equation (1) by multiplying (6) and (7) on the right and on the left by Q^{1/2}. Thus, if A is nonsingular, the minimal solution of (1) is X_{−} = Q − Y_{+}, where Y_{+} is the maximal solution of the equation

Y + AY^{−1}A^{∗} = Q.  (8)

Let r(T) be the numerical radius of an m × m matrix T, defined as

r(T) = max{|x^{∗}T x| : x ∈ C^{m}, x^{∗}x = 1}.

Concerning the spectral properties of X_{+}^{−1}A, it is proved in [28] that ρ(X_{+}^{−1}A) ≤ 1. In [15] the authors characterize the eigenvalues of X_{+}^{−1}A and give necessary and sufficient conditions so that ρ(X_{+}^{−1}A) < 1:

Theorem 2.3 The eigenvalues of X_{+}^{−1}A are precisely the eigenvalues of the matrix pencil λF − H, inside the closed unit disk, with half of the partial multiplicities for each eigenvalue on the unit circle, where

F = [ \begin{matrix}I & 0 & 0 \\ 0 & 0 & 0 \\ 0 & −I & 0 \end{matrix} ] , H = [ \begin{matrix}0 & 0 & −I \\ Q & −I & A^{∗} \\ −A & 0 & 0 \end{matrix} ] .

From the above theorem it follows that the eigenvalues of X_{+}^{−1}A are the zeros in the closed unit disk of the polynomial det(−A + λQ − λ^{2}A^{∗}).

Theorem 2.4 It holds ρ(X_{+}^{−1}A) < 1 if and only if r(Q^{−1/2}AQ^{−1/2}) < 1/2.

Moreover, in [10] the following necessary and sufficient condition is given:

Theorem 2.5 It holds ψ(λ) > 0 for all λ ∈ C if and only if r(Q^{−1/2}AQ^{−1/2}) < 1/2.

From the above theorems it follows that X_{+}^{−1}A has spectral radius strictly less than one if and only if ψ(λ) is positive definite on the unit circle.

2.2 The equation X − A^{∗}X^{−1}A = Q

Concerning equation (2), in [11] the authors prove the following interesting results:

Theorem 2.6 The rational matrix function

ψ(λ) = λA + Q − λ^{−1}A^{∗}  (9)

is nonsingular for any λ ∈ C.

Theorem 2.7 The set of solutions of (2) is non-empty, and admits a maximal element X_{+} and a minimal element X_{−}. X_{+} is the unique positive definite solution, and if A is nonsingular, X_{−} is the unique negative definite solution. Moreover, if A is nonsingular, X_{−} = Q − Y_{+}, where Y_{+} is the maximal solution of the equation

Y − AY^{−1}A^{∗} = Q.  (10)

Concerning the matrix X_{+}^{−1}A, in [11] it is proved that its spectral radius is strictly less than one, and in [15] the following characterization of the eigenvalues is given:

Theorem 2.8 The eigenvalues of X_{+}^{−1}A are precisely the eigenvalues of the matrix pencil λF − H inside the unit disk, where

F = [ \begin{matrix}I & 0 & 0 \\ 0 & 0 & 0 \\ 0 & −I & 0 \end{matrix} ] , H = [ \begin{matrix}0 & 0 & −I \\ Q & −I & −A^{∗} \\ −A & 0 & 0 \end{matrix} ] .

3 Computation of the extreme solutions of X + A^{∗}X^{−1}A = Q

In this section we describe the new algorithm for the computation of the extreme solutions of (1) and compare it with the known methods.

Throughout this section we suppose that ψ(λ) is regular and positive semidefinite for any λ ∈ C, where ψ(λ) is given in (5), so that the conditions of Theorem 2.1 are satisfied.

Let X be a solution of (1). Then, by multiplying on the right both sides of (1) by X^{−1}, we find that

−I + QX^{−1} − A^{∗}X^{−1}AX^{−1} = 0.  (11)

Thus, the matrix G = X^{−1}A solves the quadratic matrix equation

−A + QG − A^{∗}G^{2} = 0.  (12)

In addition, if Y is a solution of the matrix equation (8) then we have

−I + QY^{−1} − AY^{−1}A^{∗}Y^{−1} = 0,  (13)

and thus H = Y^{−1}A^{∗} solves the quadratic matrix equation

−A^{∗} + QH − AH^{2} = 0.  (14)

In particular, if X_{+} is the maximal solution of (1) and Y_{+} is the maximal solution of (8), then the matrix equations (12), (14) have a solution G_{+} = X_{+}^{−1}A, H_{+} = Y_{+}^{−1}A^{∗}, respectively, with spectral radius at most 1. Moreover, for Theorem 2.3, the eigenvalues of G_{+} and H_{+} are

λ(G_{+}) = {λ : g(λ) = 0, |λ| ≤ 1}

and

λ(H_{+}) = {λ : h(λ) = 0, |λ| ≤ 1},

respectively, where g(λ) = det(−A+Qλ−A^{∗}λ^{2}), h(λ) = det(−A^{∗} +Qλ−Aλ^{2}). Since g(λ) = 0 if and only if g(λ)^{∗} = 0, and since

g(λ)^{∗} = det(−A^{∗} + Qλ̄ − Aλ̄^{2}) = h(λ̄),

we conclude that λ is an eigenvalue of G_{+} if and only if λ̄ is an eigenvalue of H_{+}. In particular it holds ρ(G_{+}) = ρ(H_{+}).

The nice relation between the matrix equations (1), (8) and the quadratic matrix equations (12), (14), and the spectral properties of the solutions of the latter equations, allow us to derive a fast algorithm for the simultaneous computation of X_{+} and Y_{+}.

If A is nonsingular, the matrix Y_{+} enables us to recover the minimal solution X_{−} of (1) by means of the relation X_{−} = Q − Y_{+}.

The matrices X_{+} and Y_{+} can be simultaneously and efficiently computed by rewriting the matrix equations (11), (12), (13), (14) in terms of infinite linear systems, and by applying the cyclic reduction algorithm, according to the ideas developed in [5, 6, 7]. In fact, we observe that the following equations are satisfied:

[ \begin{matrix}Q & −A^{∗} & & 0 \\ −A & Q & −A^{∗} \\ & −A & Q & ... \\ 0 & & ... & ... \end{matrix} ] [ I ; G_{+} ; G_{+}^{2} ; ⋮ ] X_{+}^{−1} = [ I ; 0 ; 0 ; ⋮ ] ,  (15)

[ \begin{matrix}Q & −A & & 0 \\ −A^{∗} & Q & −A \\ & −A^{∗} & Q & ... \\ 0 & & ... & ... \end{matrix} ] [ I ; H_{+} ; H_{+}^{2} ; ⋮ ] Y_{+}^{−1} = [ I ; 0 ; 0 ; ⋮ ] .  (16)

By following the strategy successfully devised in [5, 6, 7] for solving nonlinear matrix equations arising in Markov chains, we apply the cyclic reduction algorithm to the above systems. This consists in performing an even-odd permutation of the block rows and columns, followed by one step of Gaussian elimination, thus generating the sequence of systems:

\begin{matrix}[ \begin{matrix}X_{n} & −A_{n}^{∗} & & 0 \\ −A_{n} & Q_{n} & −A_{n}^{∗} \\ & −A_{n} & Q_{n} & ... \\ 0 & & ... & ... \end{matrix} ] [ I ; G_{+}^{2^{n}} ; G_{+}^{2·2^{n}} ; ⋮ ] X_{+}^{−1} = [ I ; 0 ; 0 ; ⋮ ] , \\ [ \begin{matrix}Y_{n} & −A_{n} & & 0 \\ −A_{n}^{∗} & Q_{n} & −A_{n} \\ & −A_{n}^{∗} & Q_{n} & ... \\ 0 & & ... & ... \end{matrix} ] [ I ; H_{+}^{2^{n}} ; H_{+}^{2·2^{n}} ; ⋮ ] Y_{+}^{−1} = [ I ; 0 ; 0 ; ⋮ ] , n ≥ 0. \end{matrix}  (17)

The block entries of each system are defined by the following recursions:

\begin{matrix}A_{0} = A, Q_{0} = X_{0} = Y_{0} = Q, \\ A_{n+1} = A_{n}Q_{n}^{−1}A_{n}, \\ Q_{n+1} = Q_{n} − A_{n}Q_{n}^{−1}A_{n}^{∗} − A_{n}^{∗}Q_{n}^{−1}A_{n}, \\ X_{n+1} = X_{n} − A_{n}^{∗}Q_{n}^{−1}A_{n}, \\ Y_{n+1} = Y_{n} − A_{n}Q_{n}^{−1}A_{n}^{∗}, n ≥ 0. \end{matrix}  (18)

Observe that the matrices Q_{n}, X_{n} and Y_{n} are Hermitian, and thus the two block tridiagonal matrices in (17) are Hermitian.

The spectral theory of Hermitian block Toeplitz matrices [8, 24, 25, 22, 21, 27] guarantees the positive definitiveness, and thus the nonsingularity, of the blocks Q_{n}. To prove this, let us introduce the function f : (−π, π) → H^{m}, defined as f(θ) = −e^{iθ}A + Q − e^{−iθ}A^{∗}, where i is the imaginary unit, and H^{m} is the set of m × m Hermitian matrices. Let us denote by

μ_{1} = inf_{θ∈(−π,π)} λ_{min}(f(θ)),

μ_{2} = sup_{θ∈(−π,π)} λ_{max}(f(θ)),

where λ_{min}(f(θ)) (resp. λ_{max}(f(θ))) is the smallest (resp. largest) eigenvalue of f(θ).

Since ψ(λ) ≥ 0 for any λ ∈ C, it holds f(θ) ≥ 0 and μ_{1} ≥ 0. Moreover, since ψ(λ) is regular and f(θ) is a trigonometric matrix polynomial, then f(θ) is positive definite for any θ in the interval (−π, π) except at most a finite number of points.

From these properties, the following result, which guarantees the applicability of the cyclic reduction algorithm and the boundness in norm of the sequence {Q_{n}}_{n}, {X_{n}}_{n} and {Y_{n}}_{n}, can be proved:

Theorem 3.1 The matrices Q_{n}, X_{n}, Y_{n}, n ≥ 0, are positive definite, and their eigenvalues belong to the interval [μ_{1}, μ_{2}]. Moreover, it holds 0 < Q_{n+1} ≤ Q_{n}, 0 < X_{n+1} ≤ X_{n}, 0 < Y_{n+1} ≤ Y_{n}, for n ≥ 0.

Proof. For the properties of cyclic reduction [7, 4, 3], the matrices Q_{n} can be viewed as a Schur complement in T_{2^{n+1}−1}, where T_{n} is the n × n block matrix obtained by truncating the matrix in (15) at the block size n; similarly X_{n} and Y_{n} can be viewed as Schur complements in T_{2^{n}} and U_{2^{n}}, respectively, where U_{n} is the n × n block matrix obtained by truncating the matrix in (16) at the block size n. On the other hand, since T_{n} and U_{n} are Hermitian block Toeplitz matrices, their eigenvalues belong to the interval [μ_{1}, μ_{2}] (see [8, 24, 25, 22, 21, 27]); moreover, since the set where f(θ) is positive definite has strictly positive Lebesgue measure, the matrices T_{n} and U_{n} are positive definite for any n [8, 24, 25, 22, 21, 27]. In particular the Schur complements Q_{n}, X_{n}, Y_{n} are also positive definite and their eigenvalues belong to the interval [μ_{1}, μ_{2}] (see [26]). From (18) it follows that the sequences {Q_{n}}_{n}, {X_{n}}_{n} and {Y_{n}}_{n} are monotone. ⊓ ⊔

If ψ(λ) > 0 for any λ ∈ C then μ_{1} > 0, thus {Q_{n}^{−1}}_{n} is bounded in norm, and the condition number of the matrices {Q_{n}}_{n} is bounded. In this case the sequences {X_{n}}_{n} and {Y_{n}}_{n} quadratically converge to X_{+} and Y_{+}, respectively:

Theorem 3.2 If ψ(λ) > 0 for any λ ∈ C, then for any ϵ > 0 and for any matrix norm || · || it holds ||I − X_{n}X_{+}^{−1}|| = O ((σ + ϵ)^{2·2^{n}}), ||I − Y_{n}Y_{+}^{−1}|| = O ((σ + ϵ)^{2·2^{n}}), where σ = ρ(X_{+}^{−1}A) < 1. Moreover, it holds ||A_{n}|| = O ((σ + ϵ)^{2^{n}}).

Proof. From (17) we obtain that, for any n, the following relations hold:

−I + X_{n}X_{+}^{−1} − A_{n}^{∗}G_{+}^{2^{n}}X_{+}^{−1} = 0,  (19)
−A_{n} + Q_{n}G_{+}^{2^{n}} − A_{n}^{∗}G_{+}^{2·2^{n}} = 0,  (20)
−I + Y_{n}Y_{+}^{−1} − A_{n}H_{+}^{2^{n}}Y_{+}^{−1} = 0,  (21)
−A_{n}^{∗} + Q_{n}H_{+}^{2^{n}} − A_{n}H_{+}^{2·2^{n}} = 0.  (22)

From (20) and (22) we have

A_{n}(I − H_{+}^{2·2^{n}}G_{+}^{2·2^{n}}) = Q_{n}(I − H_{+}^{2^{n}}G_{+}^{2^{n}})G_{+}^{2^{n}}.  (23)

Since ρ(G_{+}) = ρ(H_{+}) = σ < 1 for Theorems 2.4 and 2.5, for any ϵ > 0 there exists an matrix norm || · || such that ||G_{+}|| ≤ σ + ϵ. Thus, since {Q_{n}}_{n} is bounded in norm, from (23) it follows that ||A_{n}|| = O ((σ + ϵ)^{2^{n}}). From (19) and (21) we obtain that ||I − X_{n}X_{+}^{−1}|| = O ((σ + ϵ)^{2·2^{n}}). Similarly, we prove that ||I − Y_{n}Y_{+}^{−1}|| = O ((σ + ϵ)^{2·2^{n}}). ⊓ ⊔

The algorithm for the computation of X_{+} and X_{−} is synthesized by the following scheme:

Algorithm 3.1 (Cyclic reduction for (1))

1 Set X_{0} = Y_{0} = Q_{0} = Q, A_{0} = A.

2 For n = 0, 1, … compute

A_{n+1} = A_{n}Q_{n}^{−1}A_{n},
Q_{n+1} = Q_{n} − A_{n}Q_{n}^{−1}A_{n}^{∗} − A_{n}^{∗}Q_{n}^{−1}A_{n},
X_{n+1} = X_{n} − A_{n}^{∗}Q_{n}^{−1}A_{n},
Y_{n+1} = Y_{n} − A_{n}Q_{n}^{−1}A_{n}^{∗},

until ||X_{n+1} − X_{n}||_{∞} < ϵ, ||Y_{n+1} − Y_{n}||_{∞} < ϵ, for a fixed error bound

ϵ > 0.

3 X_{n+1} provides an approximation to X_{+}, and, if A is nonsingular, Q − Y_{n+1} provides an approximation to X_{−}.

If we are interested only in the computation of the maximal solution X_{+} it is sufficient to skip the computation of the sequence {Y_{n}}_{n}.

One step of the above algorithm requires the solution of two Hermitian linear systems (i.e., the computation of A_{n}Q_{n}^{−1} and A_{n}^{∗}Q_{n}^{−1}), with bounded condition number, and the computation of three matrix products.

If the hypothesis of Theorem 3.2 is not satisfied, then ρ(X_{+}^{−1}A) = 1. In this case {X_{n}}_{n}, {Y_{n}}_{n}, {Q_{n}}_{n} are still Hermitian positive definite and bounded in norm. In general, the sequence {Q_{n}^{−1}}_{n} may be unbounded. However, if the sequence {A_{n}}_{n} converges to zero and the sequence {(X_{+}^{−1}A)^{2^{n}}}_{n} is bounded, the sequence {X_{n}}_{n} still converges to X_{+}.

Our algorithm can be compared with classical algorithms, which are based on fixed point iterations or Newton’s method. The basic fixed point iteration consists in the following:

Algorithm 3.2 (Fixed point iteration for (1))

1 Set X_{0} = Q.

2 For n = 0, 1, 2, …, compute X_{n+1} = Q − A^{∗}X_{n}^{−1}A until ||X_{n+1} − X_{n}||_{∞} < ϵ.

3 X_{n+1} provides an approximation to X_{+}.

In [15] it is proved that \{X_{n}\}_{n} converges to X_{+} with lim sup_{n→∞} ^{n} \sqrt{||X_{n} − X_{+}||} ≤ ρ(X_{+}^{−1}A)^{2}. Thus the sequence generated by Algorithm 3.2 converges linearly when ρ(X_{+}^{−1}A) < 1. Each step of the above algorithm requires the solution of an Hermitian linear system and a matrix product; hence the cost of the cyclic reduction algorithm is roughly twice that of the fixed point iteration. On the other hand, if we wish to compute also the minimal solution X_{−} with fixed point iterations, we must apply from scratch the algorithm to the matrix equation (8), or apply an “ad hoc” fixed point iteration described in [10]; thus the computational cost per step is the same as cyclic reduction.

To avoid possible numerical instability problems due to matrix inversions, the following inversion free variant is proposed in [15]:

Algorithm 3.3 (Inversion free fixed point iteration for (1))

1 Set X_{0} = Q, Y_{0} = I/||Q||_{∞}.

2 For n = 0, 1, 2, …, compute

Y_{n+1} = Y_{n}(2I − X_{n}Y_{n}),
X_{n+1} = Q − A^{∗}Y_{n+1}A,
until ||X_{n+1} − X_{n}||_{∞} < ϵ.

3 X_{n+1} provides an approximation to X_{+}.

The latter algorithm requires four matrix products at each step; in [15] it is shown that its convergence rate is roughly the same as that of Algorithm 3.2.

Finally, Newton’s method consists in the following:

Algorithm 3.4 (Newton’s method for (1))

1 Set X_{0} = Q.

2 For n = 1, 2, …, compute L_{n} = X_{n−1}^{−1}A and solve

X_{n} − L_{n}^{∗}X_{n}L_{n} = Q − 2L_{n}^{∗}A,  (24)
until ||X_{n} − X_{n−1}||_{∞} < ϵ.

3 X_{n} provides an approximation to X_{+}.

Each step of Newton’s method requires the solution of the Hermitian system X_{n−1}L_{n} = A, one matrix product and the solution of the Stein equation (24); thus the computational cost of each step is much larger than the cost of fixed point iterations and cyclic reduction. If equation (24) is solved with the algorithm described in [12] the computational time for each iteration is roughly 10-15 times that for Algorithm 3.2 (see [15]). Moreover, if X_{+}^{−1}A has eigenvalues very close to the unit circle, equation (24) can be nearly singular [15]. In [15] it is proved that the convergence of {X_{n}}_{n} to X_{+} is quadratic if ρ(X_{+}^{−1}A) < 1, either quadratic or linear with rate 1/2 if ρ(X_{+}^{−1}A) = 1 and all the eigenvalues of X_{+}^{−1}A on the unit circle are semi-simple.

Hence Newton’s method presents the same nice features of cyclic reduction in terms of convergence, but it is much more expensive in terms of computational cost, and more difficult to be implemented.

4 Computation of the extreme solutions of X − A^{∗}X^{−1}A = Q

For the computation of the maximal and minimal solutions X_{+} and X_{−}, respectively, of X − A^{∗}X^{−1}A = Q we can apply a technique similar to the one used in the previous section. Specifically, we introduce the maximal solution Y_{+} of (10), and observe that it holds:

−I + QX_{+}^{−1} + A^{∗}X_{+}^{−1}AX_{+}^{−1} = 0,
−I + QY_{+}^{−1} + AY_{+}^{−1}A^{∗}Y_{+}^{−1} = 0.

Thus, by setting G_{+} = X_{+}^{−1}A and H_{+} = Y_{+}^{−1}A^{∗}, the following linear systems are verified:

\begin{matrix}[ \begin{matrix}Q & A^{∗} & & 0 \\ −A & Q & A^{∗} \\ & −A & Q & ... \\ 0 & & ... & ... \end{matrix} ] [ I ; G_{+} ; G_{+}^{2} ; ⋮ ] X_{+}^{−1} = [ I ; 0 ; 0 ; ⋮ ] , \\ [ \begin{matrix}Q & A & & 0 \\ −A^{∗} & Q & A \\ & −A^{∗} & Q & ... \\ 0 & & ... & ... \end{matrix} ] [ I ; H_{+} ; H_{+}^{2} ; ⋮ ] Y_{+}^{−1} = [ I ; 0 ; 0 ; ⋮ ] . \end{matrix}  (25)

The infinite matrices in the above systems are block Toeplitz, but are not Hermitian. However, if we apply one step of cyclic reduction, we obtain the following systems:

\begin{matrix}[ \begin{matrix}X_{1} & −A_{1}^{∗} & & 0 \\ −A_{1} & Q_{1} & −A_{1}^{∗} \\ & −A_{1} & Q_{1} & ... \\ 0 & & ... & ... \end{matrix} ] [ I ; G_{+}^{2} ; G_{+}^{4} ; ⋮ ] X_{+}^{−1} = [ I ; 0 ; 0 ; ⋮ ] , \\ [ \begin{matrix}Y_{1} & −A_{1} & & 0 \\ −A_{1}^{∗} & Q_{1} & −A_{1} \\ & −A_{1}^{∗} & Q_{1} & ... \\ 0 & & ... & ... \end{matrix} ] [ I ; H_{+}^{2} ; H_{+}^{4} ; ⋮ ] Y_{+}^{−1} = [ I ; 0 ; 0 ; ⋮ ] , \end{matrix}  (26)

where

\begin{matrix}A_{1} = AQ^{−1}A, \\ Q_{1} = Q + AQ^{−1}A^{∗} + A^{∗}Q^{−1}A, \\ X_{1} = Q + A^{∗}Q^{−1}A, Y_{1} = Q + AQ^{−1}A^{∗}. \end{matrix}  (27)

Thus, after one step of cyclic reduction, we obtain two systems with Hermitian matrices, with the structure of (17), where the diagonal blocks are positive definite matrices. Hence, one step of cyclic reduction leads to an Hermitian system, like in the case of equation (1). Moreover, observe that the function

ψ_{1}(λ) = λA_{1} + Q_{1} + λ^{−1}A_{1}^{∗},

is such that ψ_{1}(λ) > 0 for any λ ∈ C, since ψ_{1}(λ) = Q+(A+λ̄A^{∗})Q^{−1}(A^{∗}+λA). If we apply cyclic reduction to systems (26), we generate the sequence (17), for n ≥ 1, where

\begin{matrix}A_{n+1} = A_{n}Q_{n}^{−1}A_{n}, \\ Q_{n+1} = Q_{n} − A_{n}Q_{n}^{−1}A_{n}^{∗} − A_{n}^{∗}Q_{n}^{−1}A_{n}, \\ X_{n+1} = X_{n} − A_{n}^{∗}Q_{n}^{−1}A_{n}, \\ Y_{n+1} = Y_{n} − A_{n}Q_{n}^{−1}A_{n}^{∗}, n ≥ 1, \end{matrix}  (28)

and A_{1}, Q_{1}, X_{1} and Y_{1} are defined in (27).

Without any assumption, the following convergence result can be proved:

Theorem 4.1 For the sequences of matrices {X_{n}}_{n}, {Y_{n}}_{n}, {Q_{n}}_{n}, defined in (28), it holds:

1. 0 < X_{n+1} ≤ X_{n}, 0 < Y_{n+1} ≤ Y_{n}, 0 < Q_{n+1} ≤ Q_{n}, n = 1, 2, …;

2. {Q_{n}}_{n} and {Q_{n}^{−1}}_{n} are bounded in norm.

Moreover, for any ϵ > 0 and for any matrix norm || · || it holds

||I − X_{n}X_{+}^{−1}|| = O ((σ + ϵ)^{2·2^{n}}) ,
||I − Y_{n}Y_{+}^{−1}|| = O ((σ + ϵ)^{2·2^{n}}) ,
||A_{n}|| = O ((σ + ϵ)^{2^{n}}) ,

where σ = ρ(X_{+}^{−1}A) < 1.

Proof. The matrices Q_{n} and A_{n}, n ≥ 2, can be seen as the matrices generated by relations (18), starting from Q_{1} and A_{1} defined in (27). Since ψ_{1}(λ) > 0 for any λ ∈ C, from Theorem 3.1 it follows that Q_{n} is positive definite for any n, and both sequences {Q_{n}}_{n} and {Q_{n}^{−1}}_{n} are bounded in norm; the monotonicity of {Q_{n}}_{n} trivially holds. The convergence properties of the sequences {X_{n}}_{n}, {Y_{n}}_{n} and {A_{n}}_{n} follow exactly as in the proof of Theorem 3.2. Since by construction X_{n} ≥ X_{n+1} (resp. Y_{n} ≥ Y_{n+1}) for any n ≥ 1, and since {X_{n}}_{n} (resp. {Y_{n}}_{n}) converges to X_{+} (resp. Y_{+}), it follows that X_{n} > 0 (resp. Y_{n} > 0) for any n. ⊓ ⊔

From the above theorem, the quadratic convergence is always guaranteed, and the condition number of Q_{n} is always bounded.

If A is nonsingular, the matrix X_{−} can be recovered from Y_{+} by means of the relation X_{−} = Q − Y_{+} (see [11]).

The resulting algorithm is resumed below:

Algorithm 4.1 (Cyclic reduction for (2))

1 Set

A_{1} = AQ^{−1}A, Q_{1} = Q + AQ^{−1}A^{∗} + A^{∗}Q^{−1}A,
X_{1} = Q + A^{∗}Q^{−1}A, Y_{1} = Q + AQ^{−1}A^{∗}.

2 For n = 1, 2, … compute

A_{n+1} = A_{n}Q_{n}^{−1}A_{n},
Q_{n+1} = Q_{n} − A_{n}Q_{n}^{−1}A_{n}^{∗} − A_{n}^{∗}Q_{n}^{−1}A_{n},
X_{n+1} = X_{n} − A_{n}^{∗}Q_{n}^{−1}A_{n},
Y_{n+1} = Y_{n} − A_{n}Q_{n}^{−1}A_{n}^{∗},

until ||X_{n+1} − X_{n}|| < ϵ, ||Y_{n+1} − Y_{n}|| < ϵ.

3 X_{n+1} provides an approximation to X_{+}, and, if A is nonsingular, Q − Y_{n+1} provides an approximation to X_{−}.

If we are interested only in X_{+}, it is sufficient to skip the computation of the sequence {Y_{n}}_{n}.

The above algorithm has the same nice features as Algorithm 3.1, in terms of computational cost and convergence properties.

We may compare our algorithm with the ones customarily used in the literature, which we recall below:

Algorithm 4.2 (Fixed point iteration for (2))

1 Set X_{0} = Q.

2 For n = 0, 1, 2, …, compute X_{n+1} = Q+A^{∗}X_{n}^{−1}A, until ||X_{n+1} −X_{n}||_{∞} < ϵ.

3 X_{n+1} provides an approximation to X_{+}.

Algorithm 4.3 (Newton’s method for (2))

1 Set X_{0} sufficiently close to X_{+}.

2 For n = 1, 2, …, compute L_{n} = X_{n−1}^{−1}A and solve

X_{n} + L_{n}^{∗}X_{n}L_{n} = Q + 2L_{n}^{∗}A,  (29)

3 X_{n} provides an approximation to X_{+}.

The above algorithms present the same characteristics as Algorithms 3.2 and 3.4, in terms of computational cost and convergence rate. Concerning Newton’s method, the initial guess X_{0} must be chosen close to the solution X_{+} in order to guarantee the convergence [15]; in practice, the starting point X_{0} is chosen as the approximation provided after a few steps of fixed point iterations (Algorithm 4.2).

The cyclic reduction algorithm has guaranteed quadratic convergence, low computational cost per iteration, and good numerical stability properties; moreover, unlike Newton’s method, it does not need an initial guess X_{0} sufficiently close to the solution X_{+}.

5 Numerical results

Concerning (1), we have compared the proposed Algorithm 3.1 with the fixed point iterations (Algorithms 3.2 and 3.3), and with Newton’s method (Algorithm 3.4), for the computation of X_{+}; for equation (2) we have compared Algorithm 4.1 with the fixed point iteration (Algorithm 4.2).

We implemented the algorithms in Fortran 90, and run the programs on a Pentium II. For the solution of systems (24) and (29) we have used the Fortran implementation of [12].

We considered the following problems:

Example 5.1 This is an example concerning (1), presented in [15]. The matrices A and Q are

A = [ \begin{matrix}0.37 & 0.13 & 0.12 \\ −0.30 & 0.34 & 0.12 \\ 0.11 & −0.17 & 0.29 \end{matrix} ] , Q = [ \begin{matrix}1.20 & −0.30 & 0.10 \\ −0.30 & 2.10 & 0.20 \\ 0.10 & 0.20 & 0.65 \end{matrix} ] .

For the stopping condition of all the algorithms we have chosen ϵ = 10^{−11}. To reach the required accuracy, the cyclic reduction (Algorithm 3.1) needs 7 iterations, fixed point iteration (Algorithm 3.2) 49 iterations, the inversion free variant (Algorithm 3.3) 52 iterations, and Newton’s method (Algorithm 3.4) 6 iterations. The CPU time needed by all the algorithms is negligible, since it is less than 0.01 seconds.

Example 5.2 We consider again equation (1), with Q = I and A symmetric, whose entries are defined by the following scheme:

1. Fix a real 0 ≤ α < 1/2.

2. For i = 1, …, m:

(a) for j = i, …, m, set a_{i,j} = i^{2} + j;

(b) compute s_{1} = ∑_{j=1}^{i−1} a_{i,j}, s_{2} = ∑_{j=i}^{m} a_{i,j};

(c) for j = i, …, m, set

a_{i,j} = a_{i,j} \frac{(1/2 − α − s_{1})}{s_{2}}, a_{j,i} = a_{i,j}.

The matrix A obtained in this way is symmetric, nonnegative and such that Ae = (1/2−α)e, where e is the vector having all the entries equal to 1. Thus the spectral radius of A is 1/2 − α, and r(Q^{−1/2}AQ^{−1/2}) = r(A) = ρ(A) = 1/2 − α. If α = 0 it holds r(A) = 1/2 and ρ(X_{+}^{−1}A) = 1; thus, the hypothesis of Theorem 3.2 are not satisfied. Small values of α provide very critical test problems.

For this example, since Q = I and A is normal with ||A||_{2} ≤ 1/2, the maximal solution is given by the matrix [29]

X_{+} = \frac{1}{2} (I + (I − 4A^{∗}A)^{1/2}) .

For the stopping condition of all the tested algorithms we have chosen ϵ = 10^{−11} when α > 0, and ϵ = 10^{−8} when α = 0. In fact, in the singular case when α = 0 a smaller value of ϵ does not provide a significantly different relative error of the computed approximation.

Table 1 reports, for different values of the size m and of the parameter α, the number of iterations needed to satisfy the stopping conditions and the relative error ||X̃_{+} − X_{+}||_{∞}/||X_{+}||_{∞}, where X̃_{+} is the approximation provided by Algorithms 3.1, 3.2, 3.3 and 3.4 (denoted by CR, FPI, IF-FPI and NM, respectively). Table 2 reports the CPU time (in seconds) needed by the same algorithms.

It is worth noting that Algorithm 3.1 shows a good convergence also in the case α = 0, i.e. ρ(X_{+}^{−1}A) = 1, when the hypothesis of Theorem 3.2 are not satisfied.

Example 5.3 This example is taken from [15] and concerns equation (2). Here

A = [ \begin{matrix}50 & 20 \\ 10 & 60 \end{matrix} ] , Q = [ \begin{matrix}3 & 2 \\ 2 & 4 \end{matrix} ] .

For the stopping condition of all the algorithms we have chosen ϵ = 10^{−11}. The cyclic reduction (Algorithm 4.1) needs 9 iterations to reach the required accuracy, while Algorithm 4.2 needs 426 iterations. The residual error of the approximation provided by Algorithm 4.1 is 2.3 · 10^{−15}, while residual error of the approximation provided by Algorithm 4.2 is 9.1 · 10^{−11}.

CR  FPI  IF-FPI  NM
α  Iter. Err.  Iter.  Err.  Iter.  Err.  Iter. Err.
m = 20
0.4  4 1.4E-15  6  1.2E-14  6  2.3E-14  3 4.0E-15
0.2  5 1.4E-15 11 2.6E-12  11 5.2E-12  4 3.7E-15
0.1  6 1.5E-15 16 2.9E-11  17 1.4E-11  5 3.9E-15
0.01  7 1.3E-15 50 1.7E-10  51 1.5E-10  6 3.1E-15
0.001  9 2.3E-15 141 6.9E-10 142 6.8E-10  8 3.2E-15
0.0001 10 5.7E-15 388 2.4E-09 389 2.4E-09  9 3.7E-15
0 26 5.3E-09 5872 5.8E-05 5873 5.9E-05 24 1.3E-08
m = 40
0.4  4 2.9E-15  6  1.2E-14  6  2.3E-14  3 6.9E-15
0.2  5 2.8E-15 11 2.6E-12  11 5.2E-12  4 6.6E-15
0.1  6 2.3E-15 16 2.9E-11  17 1.4E-11  5 5.9E-15
0.01  7 1.8E-15 50 1.7E-10  51 1.5E-10  6 5.3E-15
0.001  9 1.9E-15 141 6.8E-10 142 6.7E-10  8 5.3E-15
0.0001 10 4.5E-15 388 2.4E-09 389 2.4E-09  9 1.2E-14
0 26 5.0E-09 5822 5.8E-05 5823 5.8E-05 25 9.5E-09
m = 80
0.4  4 5.1E-15  6 1.31E-14  6 1.45E-14  5 4.3E-15
0.2  5 4.4E-15 11 2.6E-12  11 5.2E-12  4 6.6E-15
0.1  6 4.4E-15 16 2.9E-11  17 1.4E-11  5 5.9E-15
0.01  7 3.7E-15 50 1.7E-10  51 1.5E-10  6 5.3E-15
0.001  9 4.9E-15 141 6.8E-10 142 6.7E-10  8 5.3E-15
0.0001 10 5.4E-15 388 2.4E-09 389 2.4E-09  9 1.2E-14
0 26 2.2E-09 5798 5.8E-05 5799 5.8E-05 25 2.4E-09

Table 1: Example 5.2: number of iterations and errors

α  CR FPI IF-FPI NM
m = 20
0.4  0.01 0.01  0.01  0.04
0.2  0.01 0.01  0.01  0.03
0.1  0.01 0.01  0.01  0.05
0.01  0.01 0.03  0.05  0.06
0.001 0.01 0.09  0.10  0.07
0.0001 0.01 0.25  0.25  0.09
0  0.01 3.50  3.59  0.21
m = 40
0.4  0.04 0.02  0.04  0.19
0.2  0.05 0.04  0.05  0.26
0.1  0.05 0.06  0.09  0.32
0.01  0.05 0.20  0.27  0.38
0.001 0.08 0.54  0.75  0.51
0.0001 0.08 1.45  2.03  0.57
0  0.23 20.00 30.18 1.50
m = 80
0.4  0.24 0.15  0.20  1.51
0.2  0.29 0.30  0.40  2.03
0.1  0.34 0.42  0.64  2.53
0.01  0.41 1.35  1.86  3.00
0.001 0.53 3.70  5.18  4.03
0.0001 0.61 10.38 14.14 4.64
0  1.41 143.26 208.34 11.92

Table 2: Example 5.2: CPU time

Acknowledgments. The author wishes to thank an anonymous referee, whose remarks have improved the whole presentation of the paper.

References

[1] N. Akar and K. Sohraby. An invariant subspace approach in M/G/1 and G/M/1 type Markov chains. Commun. Statist. Stochastic Models, 13:381– 416, 1997.

[2] W. N. Anderson Jr., T. D. Morley, and G. E. Trapp. Positive solutions to X = A − BX^{−1}B^{∗}. Linear Algebra Appl., 134:53–62, 1990.

[3] D. A. Bini, L. Gemignani, and B. Meini. Computations with infinite Toeplitz matrices and polynomials. Submitted for publication, 2000.

[4] D. A. Bini, L. Gemignani, and B. Meini. Factorization of analytic functions by means of Koenig’s theorem and Toeplitz computations. Numer. Math., 2000. To appear.

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

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

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

[8] H. Dym. Hermitian block Toeplitz matrices, orthogonal polynomials, reproducing kernel pontryagin spaces, interpolation and extension. Oper. Theory, Adv. Appl., 34:79–135, 1998. Orthogonal matrix-valued polynomials and applications, Pap. Semin. Oper. Theory, Tel Aviv/Isr.

[9] J. C. Engwerda. On the existence of a positive definite solution of the matrix equation X + A^{T} X^{−1}A = I. Linear Algebra Appl., 194:91–108, 1993.

[10] J. C. Engwerda, A. C. M. Ran, and A. L. Rijkeboer. Necessary and sufficient conditions for the existence of a positive definite solution of the matrix equation X + A^{∗}X^{−1}A = Q. Linear Algebra Appl., 186:255–275, 1993.

[11] A. Ferrante and B. C. Levy. Hermitian solutions of the equation X = Q + NX^{−1}N^{∗}. Linear Algebra Appl., 247:359–373, 1996.

[12] J. D. Gardiner, A. J. Laub, J. J. Amato, and C. B. Moler. Solution of the Sylvester matrix equation AXB^{T} + CXD^{T} = E. ACM Trans. Math. Software, 18:223–231, 1992.

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

[14] C.-H. Guo. Newtons’s method for discrete algebraic Riccati equations when the closed-loop matrix has eigenvalues on the unit circle. SIAM J. Matrix Anal. Appl., 20:279–294, 1999.

[15] C.-H. Guo and P. Lancaster. Iterative solution of two matrix equations. Math. Comp., 68:1589–1603, 1999.

[16] G. Latouche. A note on two matrices occurring in the solution of quasibirth-and-death processes. Commun. Statist. Stochastic Models, 3:251–257, 1987.

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

[18] G. Latouche and V. Ramaswami. Introduction to Matrix Analytic Methods in Stochastic Modeling. ASA-SIAM Series on Statistics and Applied Probability 5. SIAM, 1999.

[19] G. Latouche and G.W. Stewart. Numerical methods for M/G/1 type queues. In W. J. Stewart, editor, Computations with Markov Chains, pages 571–581. Kluwer Academic Publishers, 1995.

[20] B. Meini. New convergence results on functional iteration techniques for the numerical solution of M/G/1 type Markov chains. Numer. Math., 78:39–58, 1997.

[21] M. Miranda and P. Tilli. Block Toeplitz matrices and preconditioning. Calcolo, 33(1-2):79–86 (1998), 1996. Toeplitz matrices: structures, algorithms and applications (Cortona, 1996).

[22] M. Miranda and P. Tilli. Asymptotic spectra of Hermitian block Toeplitz matrices and preconditioning results. SIAM J. Matrix Anal. Appl., 21(3):867–881 (electronic), 2000.

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

[24] S. Serra. Asymptotic results on the spectra of block Toeplitz preconditioned matrices. SIAM J. Matrix Anal. Appl., 20(1):31–44 (electronic), 1999.

[25] S. Serra. Spectral and computational analysis of block Toeplitz matrices having nonnegative definite matrix-valued generating functions. BIT, 39(1):152–175, 1999.

[26] R. L. Smith. Some interlacing properties of the Schur complement of a Hermitian matrix. Linear Algebra Appl., 177:137–144, 1992.

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

[28] X. Zhan. Computing the extremal positive definite solutions of a matrix equation. SIAM J. Sci. Comput., 17:1167–1174, 1996.

[29] X. Zhan and J. Xie. On the matrix equation X + A^{T} X^{−1}A = I. Linear Algebra Appl., 247:337–345, 1996.