Solving quadratic matrix equations arising in random walks in the quarter plane

TL;DR

This paper develops efficient algorithms leveraging Toeplitz structure to solve quadratic matrix equations in queueing models, providing convergence analysis and numerical validation for computing minimal nonnegative solutions.

Contribution

It introduces specialized fixed point iteration algorithms with acceleration strategies for quadratic matrix equations with Toeplitz structure, including a structured perturbation analysis.

Findings

Algorithms effectively compute the minimal nonnegative solution G.

Numerical experiments demonstrate the efficiency and accuracy of the proposed methods.

Structured perturbation analysis provides insights into solution stability.

Abstract

Quadratic matrix equations of the kind $A_1X^2+A_0X+A_{-1}=X$ are encountered in the analysis of Quasi--Birth-Death stochastic processes where the solution of interest is the minimal nonnegative solution $G$. In many queueing models, described by random walks in the quarter plane, the coefficients $A_1,A_0,A_{-1}$ are infinite tridiagonal matrices with an almost Toeplitz structure. Here, we analyze some fixed point iterations, including Newton's iteration, for the computation of $G$ and introduce effective algorithms and acceleration strategies which fully exploit the Toeplitz structure of the matrix coefficients and of the current approximation. Moreover, we provide a structured perturbation analysis for the solution $G$. The results of some numerical experiments which demonstrate the effectiveness of our approach are reported.

Tables

Table 1. Table 1: Parameters defining the matrices A − 1 , A 0 , A 1 subscript 𝐴 1 subscript 𝐴 0 subscript 𝐴 1 A_{-1},A_{0},A_{1} in the 2-node Jackson network of [ 24 ]

Case	${\lambda }_1$	${\lambda }_2$	${\mu }_1$	${\mu }_2$	$p$	$q$
1	1	0	1.5	2	1	0
2	1	0	2	1.5	1	0
3	0	1	1.5	2	0	1
4	0	1	2	1.5	0	1
5	1	1	2	2	0.1	0.8
6	1	1	2	2	0.8	0.1
7	1	1	2	2	0.4	0.4
8	1	1	10	10	0.5	0.5
9	1	5	10	15	0.4	0.9
10	5	1	15	10	0.9	0.4

Table 2. Table 2: Conditioning of the matrix equation for the 2-node Jackson network of [ 24 ] : actual perturbations in the solution and in the Toeplitz part and related upper bounds.

Problem	Conditioning	${\Vert {\delta }_g\Vert }_w$	${\delta }_g$-bound	${\Vert {\Delta }_G\Vert }_{\infty }$	${\Delta }_G$-bound
1	9.0000	3.0601e-09	1.3211e-08	2.3583e-09	1.9817e-08
2∗	4.5000	1.5565e-09	4.8528e-09	1.5649e-09	4.8528e-09
3	4.5000	2.9937e-09	9.3015e-09	4.2850e-09	2.0256e-08
4	9.0000	4.7978e-09	2.0256e-08	1.8119e-08	2.5845e-09
5	7.5000	5.9099e-09	1.2456e-09	1.1106e-09	5.9334e-09
$6^{\ast }$	7.5000	4.4110e-10	8.2441e-09	3.8574e-10	8.3809e-09
7	30.0000	5.4766e-09	6.8300e-08	5.4809e-09	8.1645e-08
8	5.5000	4.9898e-10	4.0333e-09	5.9838e-10	4.9549e-09
9	5.1667	7.7413e-10	2.8017e-09	7.7413e-10	3.1621e-09
${10}^{\ast }$	5.1667	6.3154e-10	3.8820e-09	5.3199e-10	4.4169e-09

Table 3. Table 3: Two-node Jackson network for Problem 7 7 7 of Table 1 : CPU time in seconds (left) and number of steps (right) required by the three fixed point iteration to arrive at a residual error at most 5.0 ⋅ 10 − 14 ⋅ 5.0 superscript 10 14 5.0\cdot 10^{-14} starting with different values of X 0 subscript 𝑋 0 X_{0} .

$0$ $I$ $T\left(g\right)$ $T\left(g\right)+v{e}_1^T$ $F_1$ 58.6 54.7 65.7 69.1 $F_2$ 37.2 34.8 40.0 28.9 $F_3$ 59.5 56.2 73.3 52.2

$0$ $I$ $T\left(g\right)$ $T\left(g\right)+v{e}_1^T$ $F_1$ 735 654 668 472 $F_2$ 466 416 421 297 $F_3$ 242 215 217 152

Table 4. Table 4: Two-node Jackson network for λ 1 = 5 subscript 𝜆 1 5 \lambda_{1}=5 , λ 2 = 0.7 subscript 𝜆 2 0.7 \lambda_{2}=0.7 , μ 1 = 2 subscript 𝜇 1 2 \mu_{1}=2 , μ 2 = 2 subscript 𝜇 2 2 \mu_{2}=2 , p = 0.5 𝑝 0.5 p=0.5 , q = 0.5 𝑞 0.5 q=0.5 : CPU time in seconds (left) and number of steps (right) required by the three fixed point iteration to arrive at a residual error at most 5.0 ⋅ 10 − 14 ⋅ 5.0 superscript 10 14 5.0\cdot 10^{-14} starting with different values of X 0 subscript 𝑋 0 X_{0} . Cyclic reduction requires 8 steps for the overall CPU time of 70.8 seconds.

$0$ $I$ $T\left(g\right)$ $T\left(g\right)+v{e}_1^T$ $F_1$ 102.8 96.0 17.5 16.0 $F_2$ 49.0 46.6 9.8 9.8 $F_3$ 4981.0 4502.0 997.0 913.0

$0$ $I$ $T\left(g\right)$ $T\left(g\right)+v{e}_1^T$ $F_1$ 806 738 103 100 $F_2$ 310 285 47 46 $F_3$ 169 149 37 35

Table 5. Table 5: Assistance from idle server model from [ 24 ] for λ 1 = 0.01 subscript 𝜆 1 0.01 \lambda_{1}=0.01 , λ 2 = 2.9 subscript 𝜆 2 2.9 \lambda_{2}=2.9 , μ 1 = 0.03 subscript 𝜇 1 0.03 \mu_{1}=0.03 , μ 2 = 2.0 subscript 𝜇 2 2.0 \mu_{2}=2.0 : CPU time in seconds (left) and number of steps (right) required by the three fixed point iterations to arrive at a residual error at most 5.0 ⋅ 10 − 14 ⋅ 5.0 superscript 10 14 5.0\cdot 10^{-14} starting with different values of X 0 subscript 𝑋 0 X_{0} . A “ ∗ * ” denotes a number of steps greater than 1000. Cyclic reduction requires 5 steps and 17 seconds of CPU time.

$0$ $I$ $T\left(g\right)$ $T\left(g\right)+v{e}_1^T$ $F_1$ * * 204.8 191.6 $F_2$ 16.4 13.7 4.6 3.9 $F_3$ 405.1 402.6 218.3 183.6

$0$ $I$ $T\left(g\right)$ $T\left(g\right)+v{e}_1^T$ $F_1$ * * 844 782 $F_2$ 42 40 10 9 $F_3$ 26 25 9 7

Table 6. Table 6: Random walk in the quarter plane: Number of iterations and CPU time in seconds. From left to right: fixed point iterations with X 0 = T ​ ( g ) subscript 𝑋 0 𝑇 𝑔 X_{0}=T(g) , fixed point iterations with X 0 = T ​ ( g ) + v ​ e 1 T subscript 𝑋 0 𝑇 𝑔 𝑣 superscript subscript 𝑒 1 𝑇 X_{0}=T(g)+ve_{1}^{T} , cyclic reduction, Newton iteration with X 0 = 0 subscript 𝑋 0 0 X_{0}=0 , combination of cyclic reduction and Newton iteration. A “ ∗ * ” denotes more than 10000 iterations and a CPU larger than 1000 seconds. Starting the iterations with X 0 = T ​ ( g ) + v ​ e 1 T subscript 𝑋 0 𝑇 𝑔 𝑣 superscript subscript 𝑒 1 𝑇 X_{0}=T(g)+ve_{1}^{T} generates sequences converging to the stochastic solution G ^ ^ 𝐺 \widehat{G} , while starting with X 0 = 0 subscript 𝑋 0 0 X_{0}=0 or applying CR, Newton iteration and their combinations generate sequences converging to the minimal (substochastic) solution G 𝐺 G .

	$F_1$	$F_2$	$F_3$	$F_1$	$F_2$	$F_3$	CR	Newton	CR+Newton
iter	*	*	*	285	205	119	23	23	15+10
CPU	*	*	*	3.1	2.3	2.6	524	3.4	0.4+1.5

Equations

$$A_1{X}^2+A_{0}X+A_{-1}=X,$$

$$A_1+X{A}_0+X^2{A}_{-1}=X,$$

$$A_i=\left[\begin{array}{ccccc}{b}_{i,0}& b_{i,1}& & & \\ a_{i,-1}& a_{i,0}& a_{i,1}& & \\ & a_{i,-1}& a_{i,0}& a_{i,1}& \\ & & \ddots & \ddots & \ddots \end{array}\right],i=-1,0,1,$$

$$a_1(z)g(z{)}^2+a_0(z)g(z)+a_{-1}(z)=g(z),|z|=1.$$

$$a_i(z)=a_{i,-1}{z}^{-1}+a_{i,0}+a_{i,1}z$$

$$b_i(z)=b_{i,0}+b_{i,1}z,$$

$$\phi (z)=-(I-zR)W(I-z^{-1}G)$$

$$A_1=RW,A_{-1}=WG,A_0=I-W-RWG,$$

$$W=I-A_0-A_{1}G=I-A_0-R{A}_{-1},$$

$$G1\le 1.$$

$$\theta =\min \{a_{-1}(1),b_{-1}(1)\},\gamma =\max \left\{\frac{a_1(1)}{a_{-1}(1)},\frac{b_1(1)}{b_{-1}(1)}\right\}.$$

$$\Vert W^{-1}{\Vert }_{\infty }=\Vert \sum _{k=0}^{\infty }(A_0+A_{1}G{)}^{k}1{\Vert }_{\infty }\le \frac{1}{1-\Vert A_0+A_1{\Vert }_{\infty }}=\frac{1}{\theta }.$$

$$W^{-1}{A}_{1}1\le \gamma W^{-1}{A}_{-1}1=\gamma G1\le \gamma 1,$$

$$\mathcal{Q}\mathcal{T}:=\{A=T(f)+E\}$$

$$a_{-1}(z)+a_0(z)\lambda +a_1(z){\lambda }^2=\lambda .$$

$$A_{-1}1>A_{1}1,or{A}_{-1}1\ge A_{1}1>0,$$

$$A_{-1}1>A_{1}1$$

$$g(z)=\sum _{j=-n+1}^n{g}_j{z}^j+\sum _{j=-n+1}^n\sum _{k\ge 1}(z^{mk+j}{g}_{mk+j}+z^{-mk-j+1}{g}_{-mk-j+1}).$$

$$g({\omega }_m^i)=\sum _{j=-n+1}^n{g}_j{\omega }_m^{ij}+\sum _{j=-n+1}^n({\omega }_m^{ij}\sum _{k\ge 1}{g}_{mk+j}+{\omega }_m^{-i(j-1)}\sum _{k\ge 1}{g}_{-mk-j+1}).$$

$$\hat{g}(z)=\sum _{j=-n+1}^n{g}_j{z}^j+\sum _{j=-n+1}^n(z^j{\hat{g}}_j^{+}+z^{-j+1}{\hat{g}}_{-j+1}^{-}),$$

$${\hat{g}}_j^{+}=\sum _{k\ge 1}{g}_{mk+j},{\hat{g}}_{-j+1}^{-}=\sum _{k\ge 1}{g}_{-mk-j+1},j=-n+1,\dots ,n,$$

$$g^{\prime \prime }(1)-{\hat{g}}^{\prime \prime }(1)\ge 2n\left(\sum _{j<-n+1}{g}_j+\sum _{j>n}{g}_j\right),$$

$$g^{\prime \prime }(1)-{\hat{g}}^{\prime \prime }(1)=\sum _{j=-n+1}^n\sum _{k\ge 1}\left(g_{mk+j}{\alpha }_{j,k}+g_{-mk-j+1}{\alpha }_{j,k}\right),$$

$$a_1(z)g(z{)}^2+(a_0(z)-1)g(z)+a_{-1}(z)=0,$$

$$a_1^{\prime }(z)g(z{)}^2+2g(z)g^{\prime }(z)a_1(z)+(a_0(z)-1)g^{\prime }(z)+a_0^{\prime }(z)g(z)+a_{-1}^{\prime }(z)=0,$$

$$g^{\prime }(1)=\frac{a_1^{\prime }(1)g(1{)}^2+a_0^{\prime }(1)g(1)+a_{-1}^{\prime }(1)}{1-2a_1(1)g(1)-a_0(1)},g(1)=\min (1,a_{-1}(1)/a_1(1)).$$

$$a_1^{\prime \prime }(1)g(1{)}^2+$$

$$+$$

$$g^{\prime \prime }(1)=$$

$$+$$

$$F_1(X)=A_{-1}+A_{0}X+A_1{X}^2,$$

$$F_2(X)=(I-A_0{)}^{-1}(A_{-1}+A_1{X}^2),$$

$$F_3(X)=(I-A_0-A_{1}X{)}^{-1}{A}_{-1},$$

$$\sigma =1-\min (a_{-1}(1)-a_1(1),b_{-1}(1)-b_1(1))<1.$$

Full text

Keywords: Matrix equations, random walks, Markov chains, Toeplitz matrices, infinite matrices, fixed point iteration, Newton iteration.

MSC: 65F30, 15A24, 60J22, 15B05

1 Introduction

Random walks in the quarter plane describe a wide variety of two-queue models with various service policies such as nonpreemptive priority, $K$-limited service, server vacation and server setup [26]. Models of this kind concern, for instance, bi-lingual call centers [30], generalized two-node Jackson networks [27], two-demand models [10], two-stage inventory queues [11], and more.

A theoretical analysis of stability, of tail decay rates and of other asymptotic properties has been carried out by several authors, in particular in [14], [22], [23], [26], and in the book [9], in which the invariant measure and the transient behavior are investigated by means of analytic and functional tools.

A different approach is based on representing a random walk in the quarter plane as a 2-dimensional Quasi–Birth-Death (QBD) stochastic process. This latter framework, based on the matrix analytic approach of [25], allows to express the invariant probability measure, and other quantities of interest for the stochastic model, in terms of a solution of suitable quadratic matrix equations. This provides a further tool for the theoretical analysis [17], [19] and paves the way for the design of effective algorithms based on the numerical solution of quadratic matrix equations.

In fact, relying on the matrix analytic theory of [25], the problem of computing the invariant probability measure of a QBD process is reduced to computing the minimal nonnegative solution $G$ and $R$ of the two matrix equations

[TABLE]

respectively, where the coefficients $A_{-1},A_{0},A_{1}$ are nonnegative matrices such that $A_{-1}+A_{0}+A_{1}$ is row-stochastic and $X$ is the unknown. We say that a matrix $X$ is nonnegative, and we write $X\geq 0$, if its entries are nonnegative. Moreover we say that a solution $X$ of a matrix equation is minimal nonnegative if $X\geq 0$ and for any other nonnegative solution $Y$ it holds $Y-X\geq 0$. For more details in this regard, we refer the reader to the books [2], [19], and [25].

In the case where the coefficients are finite dimensional, several algorithms have been introduced to compute $G$ and $R$. They include fixed point iterations and doubling algorithms like Logarithmic Reduction and Cyclic Reduction (CR) [7], [2], [19].

In the case of 2-dimensional QBDs the coefficients $A_{-1},A_{0}$, $A_{1}$ are semi-infinite and have a special structure, more precisely,

[TABLE]

where $a_{i,j}\geq 0$, $b_{i,j}\geq 0$ and $\sum_{i,j=-1}^{1}a_{i,j}=1$, $\sum_{i=-1}^{1}\sum_{j=0}^{1}b_{i,j}=1$. These blocks belong to the class of matrices representable in the form $A=T(s)+E$ where $T(s)$ is the Toeplitz matrix associated with the symbol $s(z)=\sum_{i\in\mathbb{Z}}s_{i}z^{i}$, that is, $(T(s))_{i,j}=s_{j-i}$, and $E=(e_{i,j})$ is such that $v_{i}=\sum_{j}|e_{i,j}|$ is finite and $\lim_{i}v_{i}=0$. The matrix $T(s)$ is called Toeplitz part while $E$ is called correction. Here $s(z)$ is a function belonging to the Wiener class $\mathcal{W}=\{f(z)=\sum_{i\in\mathbb{Z}}f_{i}z^{i},~{}\|f\|_{w}:=\sum_{i\in\mathbb{Z}}|f_{i}|<\infty\}$. In particular, for the matrix in (3) it is easy to check that $A_{i}=T(a_{i})+E_{i}$ where $a_{i}(z)=\sum_{j=-1}^{1}a_{i,j}z^{j}$ and $E_{i}$ has zero entries except for the first row which is equal to $[b_{i,0}-a_{i,0},b_{i,1}-a_{i,1},0,\ldots]$. Matrices of this kind are called Quasi-Toeplitz (QT) in [3].

The case of QBD with infinite blocks has been initially investigated in [17], [18], and [20], by reducing the problem to finite size relying on truncation and augmentation of the blocks. However, this approach does not lead to reliable computational techniques since the result of the numerical computation is strongly dependent on the way the infinite matrices have been truncated. In fact, the models obtained by truncating the infinite dimensional problem may have asymptotic properties, like the decay rate, which are not consistent with the original problem [15], [29].

More recently, conditions under which the solution $G$ of (1) can be represented as the sum $G=T(g)+E_{g}$ are given in [5], so that, despite the solution $G$ has infinitely many entries, it can be represented up to any arbitrary approximation error by using a finite number of parameters. Moreover, in [3] and [4], by using the structure properties of QT matrices, the algorithm of Cyclic Reduction has been extended to the case of infinite matrices. This algorithm still keeps a fast convergence speed in terms of number of iterative steps. However, in certain cases the cost of each step becomes extremely large due to the cost of certain operations with QT matrices, like matrix inversion and the compression of the correction part. Another drawback of CR is that this iteration is not self-correcting.

In this paper, we propose and analyze some fixed point iterations which have a low cost per step and, unlike CR, are self-correcting and allow to keep separated the computation of the Toeplitz part $T(g)$ and the correction part $E_{g}$. In fact, we show that the symbol $g(z)$ defining the Toeplitz part satisfies the functional equation

[TABLE]

We use this property to design an algorithm based on evaluation and interpolation at the roots of 1 for approximating the coefficients of $g(z)$, which is extremely fast and allows for an automatic control of the number of interpolation points, according to the desired approximation error.

The correction part $E_{g}$ is obtained by simply applying fixed point iterations. We consider three iterations of the kind $X_{k+1}=F(X_{k})$, $k=0,1,\ldots$, defined by suitable functions $F_{1}$, $F_{2}$, $F_{3}$, where $F_{1}$ requires no matrix inversion, $F_{2}$ requires to compute an inverse matrix once for all, while $F_{3}$ requires one inversion per step. These iterations are well known in the case of QBD with finitely many phases [2], [19], [21], and are here extended to coefficients with the QT structure (3).

We show that, under mild assumptions, starting with $X_{0}=0$, the sequences generated by $F_{1},F_{2},F_{3}$ converge monotonically and linearly to $G$ in the infinity norm. Moreover, we prove that the rate of convergence of the sequence generated by $F_{3}$ is better than the rate of the sequence generated by $F_{2}$, which in turn is better than that generated by $F_{1}$. We prove that if $X_{0}$ is row-stochastic then all the matrices $X_{k}$ are row-stochastic and the rate of convergence of each of the three iterations is better than that obtained with $X_{0}=0$. Numerical experiments show also the evidence that for $X_{0}=T(g)$ the rate of convergence of the three sequences is even better.

Then we adapt Newton’s iteration, in the form given by [16], to the case of QT coefficients. In order to solve the Sylvester equation arising at each step of Newton’s iteration we rely on the solver introduced in [28]. Under mild conditions, we prove that, for $X_{0}=0$, convergence holds in the infinity norm, is monotonic and quadratic.

In order to evaluate an a posteriori bound on the approximation error in the computation of $G$, we also perform the analysis of the structured condition number. More specifically, we provide perturbation results related to perturbations of the Toeplitz part and of the correction part in the matrix coefficients $A_{i}$, $i=-1,0,1$, and we estimate the consequent variation of the solution $G$ and of its Toeplitz part $T(g)$.

Numerical experiments are reported which show the effectiveness of our approach and the reliability of our algorithms with respect to the algorithm CR. In particular we show that in certain cases the combination of Newton’s iteration and cyclic reduction provides a substantial acceleration of the convergence.

The paper is organized as follows: in Section 2 we recall some preliminary properties and concepts useful for the analysis of the problem; in Section 3 we present an algorithm for computing the Toeplitz part $T(g)$ of the solution; Section 4 deals with the analysis of three fixed-point iterations applied to infinite QT matrices, while Section 5 concerns the algorithmic analysis of Newton iteration; in Section 6 we carry out the analysis of the conditioning by providing some perturbation results, while in Section 7 we present and discuss some numerical experiments which show the effectiveness of our approach.

2 Preliminaries

Let $\ell^{\infty}$ be the set of sequences $x=(x_{i})_{i\in\mathbb{N}}$ such that $\|x\|_{\infty}:=\sup_{i}|x_{i}|$ is finite. Consider the set of matrices $A=(a_{i,j})$ such that the application $x\to y=Ax$, where $y_{j}=\sum_{j=1}^{\infty}a_{i,j}x_{j}$, defines a linear operator from $\ell^{\infty}$ to $\ell^{\infty}$. Denote this set by $\mathcal{L}_{\infty}$ and define the induced norm $\|A\|_{\infty}=\sup_{\|x\|_{\infty}=1}\|Ax\|_{\infty}$. It can be verified that $\|A\|_{\infty}=\sup_{i}\sum_{j=1}^{\infty}|a_{i,j}|$. Recall that $\mathcal{L}_{\infty}$ is a Banach algebra, that is, it is closed under the row-by-column product, the norm satisfies $\|AB\|_{\infty}\leq\|A\|_{\infty}\cdot\|B\|_{\infty}$ for any $A,B\in\mathcal{L}_{\infty}$, and the normed space is complete.

We introduce the following notation

[TABLE]

so that we may write $A_{i}=T(a_{i})+E_{i}$, for $i=-1,0,1$, where $E_{i}$ has null entries except for those in the first row which are equal to $[b_{i,0}-a_{i,0},b_{i,1}-a_{i,1},0,\ldots]$. We assume that the entries of $A_{i}$ are nonnegative and $(A_{-1}+A_{0}+A_{1}){\bm{1}}={\bm{1}}$, where ${\bm{1}}$ is the vector of all ones of appropriate dimension. It is known [19], [31] that under these conditions, there exist the minimal nonnegative solutions $R$ and $G$ of (1) and (2), respectively, and the Laurent matrix polynomial $\varphi(z)=z^{-1}A_{-1}+A_{0}-I+zA_{1}$ admits the factorization

[TABLE]

where

[TABLE]

Observe that if $a_{-1}(1)=0$, i.e., $A_{-1}=e_{1}w^{T}$, $w^{T}=(b_{-1,0},b_{-1,1},0,\ldots)\neq 0$, then the minimal nonnegative solution $G$ of equation (1) can be expressed in the form $G={\bm{1}}v^{T}$ where $v=\frac{1}{b_{-1}(1)}w$. Therefore, without loss of generality we may assume that $a_{-1}(1)>0$.

The following result is valid if $a_{-1}(1)>0$ and $b_{-1}(1)>0$, that is, $A_{-1}{\bm{1}}>0$.

Lemma 2.1**.**

Assume $A_{-1}{\bm{1}}>0$ and define

[TABLE]

Then the matrix $W=I-A_{0}-A_{1}G$ is invertible in $\mathcal{L}_{\infty}$, has nonnegative inverse, and $\|W^{-1}\|_{\infty}\leq\frac{1}{1-\|(A_{0}+A_{1}){\bm{1}}\|_{\infty}}=\frac{1}{\theta}$. Moreover, $\|W^{-1}RW\|_{\infty}\leq\gamma$. If $A_{-1}{\bm{1}}>A_{1}{\bm{1}}$ then $\gamma<1$.

Proof 2.2**.**

Observe that $\|A_{0}+A_{1}G\|_{\infty}=\|(A_{0}+A_{1}G){\bm{1}}\|_{\infty}\leq\|(A_{0}+A_{1}){\bm{1}}\|_{\infty}=\|(I-A_{-1}){\bm{1}}\|_{\infty}=1-\theta<1$, since $(A_{0}+A_{1}){\bm{1}}=(I-A_{-1}){\bm{1}}$, $A_{-1}{\bm{1}}>0$ and $A_{-1}{\bm{1}}=(b_{-1}(1),a_{-1}(1),a_{-1}(1),\ldots)^{T}$. Therefore $W^{-1}\in\mathcal{L}_{\infty}$ and is nonnegative, being $W^{-1}=\sum_{i=0}^{\infty}(A_{0}+A_{1}G)^{i}$ and

[TABLE]

Now we show that $\|W^{-1}RW\|_{\infty}\leq\gamma$. Since $W^{-1}RW=W^{-1}A_{1}$ by (5), it is sufficient to consider $\|W^{-1}A_{1}\|_{\infty}=\|W^{-1}A_{1}{\bm{1}}\|_{\infty}$. By definition of $\gamma$ we have $A_{1}{\bm{1}}\leq\gamma A_{-1}{\bm{1}}$, therefore

[TABLE]

*where we used the fact that $W^{-1}A_{-1}=G$. Thus we have $\|W^{-1}RW\|_{\infty}\leq\gamma$. Since $A_{1}{\bm{1}}=(b_{1}(1),a_{1}(1),a_{1}(1),\ldots)^{T}$ and $A_{-1}{\bm{1}}=(b_{-1}(1),a_{-1}(1),a_{-1}(1),\ldots)^{T}$ then $A_{-1}{\bm{1}}>A_{1}{\bm{1}}$ implies $\gamma<1$. *

Define $\mathcal{W}=\{f(z)=\sum_{i\in\mathbb{Z}}f_{i}z^{i}:\quad\|f\|_{w}:=\sum_{i\in\mathbb{Z}}|f_{i}|<\infty\}$. Consider the following class

[TABLE]

where $f(z)\in\mathcal{W}$, the matrix $E=(e_{i,j})\in\mathcal{L}_{\infty}$ is such that $\lim_{i}v_{i}=0$, where $v_{i}=\sum_{j=1}^{\infty}|e_{i,j}|$.

Observe that $A_{i}\in\mathcal{QT}$ for $i=-1,0,1$, moreover, in [5] it is shown that $\mathcal{QT}$ is an algebra with the infinity norm, and the matrices $W,G$ and $R$ in (5) belong to $\mathcal{QT}$ if $A_{-1}{\bm{1}}>A_{1}{\bm{1}}$ or if $A_{-1}{\bm{1}}\geq A_{1}{\bm{1}}>0$. More precisely we have the following result

Theorem 2.3**.**

The minimal nonnegative solution $G$ of the matrix equation (1) can be written as $G=T(g)+E_{g}$ where $g(z)=\sum_{i\in\mathbb{Z}}g_{i}z^{i}\in\mathcal{W}$ is such that $g_{i}\geq 0$ and $\|g\|_{w}=g(1)\leq 1$. Moreover, for any $z$ such that $|z|=1$, $g(z)$ is a solution of minimum modulus of the quadratic equation

[TABLE]

This solution is unique if there exists $j$ such that $a_{i,j}\neq 0$ for at least two different values of $i$. If

[TABLE]

*then $G\in\mathcal{QT}$, $G{\bm{1}}={\bm{1}}$ and $g(1)=1$. Conversely, if $G{\bm{1}}={\bm{1}}$ and $G\in\mathcal{QT}$ then $a_{-1}(1)\geq a_{1}(1)$ and $g(1)=1$. *

Observe that the condition

[TABLE]

is equivalent to $a_{-1}(1)>a_{1}(1)$ and $b_{-1}(1)>b_{1}(1)$. In the following we assume that (8) holds and that there exists $i$ such that $a_{i,j}\neq 0$ for at least two values of $j$. The latter condition is very mild.

3 Computing the symbol $g(z)$

Relying on Theorem 2.3, we provide an algorithm, based on the evaluation/interpolation at the roots of 1, for computing an approximation $\hat{g}_{i}$, $i=-n+1,\ldots,n$, to the coefficients $g_{i}$ of $g(z)$, where $n$ is such that $|\hat{g}_{i}-g_{i}|\leq\epsilon/(2n)$, for any $i$ and for a given tolerance $\epsilon>0$. In this analysis we may relax the assumption $a_{-1}(1)>a_{1}(1)$ so that the result holds in general.

Let $n>0$ be an integer and set $m=2n$. Define $\omega_{m}=\cos\frac{2\pi}{m}+\mathbf{i}\sin\frac{2\pi}{m}$ a principal $m$th root of 1, where $\mathbf{i}$ is the imaginary unit such that $\mathbf{i}^{2}=-1$. Rewrite $g(z)$ as

[TABLE]

Since $\omega_{m}^{km}=1$, from (9) we have

[TABLE]

Therefore, the Laurent polynomial defined by

[TABLE]

is such that $g(\omega_{m}^{i})=\hat{g}(\omega_{m}^{i})$, that is, it interpolates $g(z)$ at the $m$-th roots of 1.

The following lemma provides a bound to the tail of the Laurent series $g(z)$ and extends to the case of Laurent series a similar property proved in [8] valid for power series.

Lemma 3.1**.**

Let $g(z)$ be the solution of minimum modulus of equation (7). Let $\hat{g}(z)=\sum_{j=-n+1}^{n}\hat{g}_{j}z^{j}$ be the Laurent polynomial interpolating $g(z)$ at the $m$-th roots of 1, i.e., such that $g(\omega_{m}^{i})=\hat{g}(\omega_{m}^{i})$, $i=-n+1,\ldots,n$, where $m=2n$. If $g^{\prime\prime}(x)\in\mathcal{W}$, then $g^{\prime\prime}(1)\geq 0$ and

[TABLE]

*moreover $0\leq\hat{g}_{j}-g_{j}\leq\frac{1}{2n}(g^{\prime\prime}(1)-\hat{g}^{\prime\prime}(1))$, for $j=-n+1,\ldots,n$. *

Proof 3.2**.**

Since the coefficients $g_{i}$ are nonnegative then also $g^{\prime\prime}(z)$ has nonnegative coefficients, moreover, since $g^{\prime\prime}(z)\in\mathcal{W}$, then the series $g^{\prime\prime}(1)$ is absolutely convergent and $g^{\prime\prime}(1)\geq 0$. Thus, from the representation (9), in view of (10), we deduce that

[TABLE]

*where $\alpha_{j,k}=(mk+j)(mk+j-1)-j(j-1)$. The inequality (11) follows from the nonnegativity of the coefficients and from the property $\alpha_{j,k}\geq m$, valid for $k\geq 1,~{}j=-n+1,\ldots,n$ which can be verified by a direct inspection. The bound on $\hat{g}_{j}-g_{j}$ follows from (11) since $\hat{g}_{j}=g_{j}+\hat{g}_{j}^{+}+\hat{g}_{j}^{-}$ and $\hat{g}^{+}_{j}+\hat{g}_{j}^{-}\leq\sum_{j<-n+1}g_{j}+\sum_{j>n}g_{j}$ in view of (10). *

Observe that $\hat{g}^{\prime\prime}(1)$ is computable once the coefficients of the polynomial $\hat{g}(z)$ have been computed. Moreover, the value of $g^{\prime\prime}(1)$ is computable even though $g(z)$ is not known. In fact, by taking the second derivative in the equation obtained by replacing $\lambda$ with $g(z)$ in (7), i.e.,

[TABLE]

for $z=1$, the value of $g^{\prime\prime}(1)$ can be easily expressed in terms of $a_{i}(1)$, $a_{i}^{\prime}(1)$ and $a_{i}^{\prime\prime}(1)$. More precisely, by taking the first derivative we obtain

[TABLE]

which yields

[TABLE]

By taking the second derivative for $z=1$, we get

[TABLE]

which yields

[TABLE]

Lemma 3.1 provides an a posteriori bound to the error in the approximation of the Laurent series $g(z)$ together with a stop condition for the following evaluation interpolation algorithm for computing the coefficients of $g(z)$.

Observe that the error bound converges to zero at least as $O(1/n)$. If the function $g(z)$ is analytic in a neighborhood of the unit circle, then its coefficients decay exponentially to zero [12] so that also the the bound on the error converges exponentially to zero. It is also interesting to observe that, for $n\to\infty$, the convergence of the coefficients of $\hat{g}(z)$ to the corresponding coefficients of $g(z)$ is monotonic.

Finally observe that the overall computational cost of this algorithm is $O(n\log n)$ arithmetic operations. In order to complete the computation of $G$ it remains to approximate the correction $E_{g}$. In view of the fact that $E_{g}$ has entries $e^{(g)}_{i,j}$ such that $v_{i}=\sum_{j}|e^{(g)}_{i,j}|$ is finite and $\lim_{i}v_{i}=0$, we can approximate $E_{g}$ with a finite number of parameter within an error bound $\epsilon$. This computation is performed by means of functional iteration and is analyzed in the next section.

4 Fixed point iterations

In this section we analyze the convergence of sequences generated by a functional iteration of the kind $X_{k+1}=F(X_{k})$, $k=0,1,\ldots$, where $F(X)$ is a matrix function such that $G=F(G)$ where $G$ is the minimal nonnegative solution of (1). More precisely, we will consider the following cases

[TABLE]

while Newton iteration is considered in the next section.

In the case where $A_{-1},A_{0},A_{1}$ are finite matrices, the convergence analysis of the sequences generated by the functions (14) has been performed in [21]. Here, we extend the results of [21] to the case of matrices of infinite size belonging to $\mathcal{L}_{\infty}$. We need the following

Lemma 4.1**.**

Let $A_{-1}{\bm{1}}>A_{1}{\bm{1}}$, and

[TABLE]

Let $H_{1}=A_{0}+A_{1}+A_{1}G$. Then $\|H_{1}\|_{\infty}\leq\sigma$, so that $\|A_{0}\|_{\infty}\leq\sigma$, $\|A_{0}+A_{1}G\|_{\infty}\leq\sigma$. Therefore $I-A_{0}$ and $I-A_{0}-A_{1}G$ are invertible and, for $H_{2}=(I-A_{0})^{-1}(A_{1}+A_{1}G)$, $H_{3}=(I-A_{0}-A_{1}G)^{-1}A_{1}$ we have

[TABLE]

Proof 4.2**.**

We have $\|A_{0}+A_{1}+A_{1}G\|_{\infty}=\|(A_{0}+A_{1}+A_{1}G){\bm{1}}\|_{\infty}$. Moreover, since by Theorem 2.3 we have $G{\bm{1}}={\bm{1}}$ and ${\bm{1}}=(A_{-1}+A_{0}+A_{1}){\bm{1}}$, then $(A_{0}+A_{1}+A_{1}G){\bm{1}}=(A_{0}+A_{1}+A_{1}){\bm{1}}={\bm{1}}-(A_{-1}-A_{1}){\bm{1}}\leq\sigma{\bm{1}}$, by definition of $\sigma$. This implies $\|A_{0}+A_{1}+A_{1}G\|_{\infty}\leq\sigma$. Since $A_{0},A_{1},G$ are nonnegative then $\|A_{0}\|_{\infty}\leq\|A_{0}+A_{1}+A_{1}G\|_{\infty}\leq\sigma$ and $\|A_{0}+A_{1}G\|_{\infty}\leq\sigma$. The matrices $I-A_{0}$ and $I-A_{0}-A_{1}G$ are invertible in $\mathcal{L}_{\infty}$ since, in general, if $I-B\in\mathcal{L}_{\infty}$ is such that $\|B\|_{\infty}<1$ then the series $\sum_{i=0}^{\infty}B^{i}$ has norm bounded by $1/(1-\|B\|_{\infty})$ and coincides with $(I-B)^{-1}$. Concerning $H_{2}$ we have $\|H_{2}\|_{\infty}=\|H_{2}{\bm{1}}\|_{\infty}$. Moreover,

[TABLE]

where we used the properties $(A_{-1}-A_{1}G){\bm{1}}>0$ and $(I-A_{0})^{-1}\geq I$. Concerning $H_{3}$, since $A_{1}{\bm{1}}=(I-A_{0}-A_{-1}){\bm{1}}=(I-A_{0}-A_{1}G-(A_{-1}-A_{1}G)){\bm{1}}$ we have

[TABLE]

*where we used the fact that $G{\bm{1}}={\bm{1}}$, $(A_{-1}-A_{1}G){\bm{1}}>0$ and $(I-A_{0}-A_{1}G)^{-1}\geq(I-A_{0})^{-1}$. *

Observe that, from Lemma 2.1 and from (5), it follows that $H_{3}=W^{-1}RW$ and $\|H_{3}\|_{\infty}=\|W^{-1}RW\|_{\infty}\leq\gamma$ where $\gamma<1$ is defined in (6). This provides a different bound on the norm of $H_{3}$. Therefore we have

[TABLE]

with $\gamma$ and $\sigma$ defined in (6) and (15), respectively.

We are ready to prove the following result which shows that the three sequences generated by (14) starting with $X_{0}=0$ monotonically converge to $G$, convergence holds in the infinity norm and is linear. Moreover, the convergence of the third iteration is faster than that of the second one, while the convergence of the second iteration is faster than that of the first one.

Theorem 4.3**.**

Assume that $A_{-1}{\bm{1}}>A_{1}{\bm{1}}$. For $i\in\{1,2,3\}$ define $X_{k+1}^{(i)}=F_{i}(X_{k}^{(i)})$, $k=0,1,2,\ldots$, where $X_{0}^{(i)}=0$ and $F_{i}(X)$ are given in (14). Then

the three sequences $\{X_{k}^{(i)}\}$ are well defined, 2. 2.

$0\leq X_{k}^{(i)}\leq X_{k+1}^{(i)}\leq G$, 3. 3.

for the error $\mathcal{E}_{k}^{(i)}=G-X_{k}^{(i)}$ we have $\|\mathcal{E}_{k+1}^{(i)}\|_{\infty}\leq\|H_{i}\|_{\infty}\|\mathcal{E}_{k}^{(i)}\|_{\infty}$, where $H_{i}$, $i=1,2,3$ are the matrices defined in Lemma 4.1, so that $\lim_{k\to\infty}\|\mathcal{E}_{k}^{(i)}\|_{\infty}=0$.

Proof 4.4**.**

The first iteration is clearly well defined. The second is well defined since, according to Lemma 4.1, the matrix $I-A_{0}$ is invertible in $\mathcal{L}_{\infty}$. The third iteration is well defined as long as the matrix $I-A_{0}-A_{1}X_{k}$ is invertible. On the other hand if $0\leq X_{k}\leq G$, the latter matrix is invertible in view of Lemma 4.1 since $\|A_{0}+A_{1}X_{k}\|_{\infty}\leq\|A_{0}+A_{1}G\|_{\infty}$. In order to prove that $0\leq X_{k}^{(i)}\leq X_{k+1}^{(i)}\leq G$ we use an induction argument. We prove it for the first iteration, i.e., for $i=1$, the same technique can be used for the other iterations. For notational simplicity we omit the superscript and write $X_{k}$ in place of $X_{k}^{(1)}$. Since for $X_{0}=0$ we have $X_{1}=A_{-1}$, so that $0\leq X_{0}\leq X_{1}$ and $G-X_{1}=G-A_{-1}=(A_{0}+A_{1}G)G\geq 0$. For the inductive step, assume that $0\leq X_{k-1}\leq X_{k}\leq G$. We first show that $0\leq X_{k}\leq X_{k+1}$. From $X_{k+1}=A_{1}X_{k}^{2}+A_{0}X_{k}+A_{-1}$ and from the property $0\leq X_{k-1}\leq X_{k}$ we get

[TABLE]

Now consider

[TABLE]

Since $G-X_{k}\geq 0$ then also $G-X_{k+1}\geq 0$.

Concerning the norm bounds to $\mathcal{E}_{k}$, for $\mathcal{E}_{k}=\mathcal{E}_{k}^{(1)}$ for the sequence defined by $F_{1}$, from (17) we obtain

[TABLE]

Since $\mathcal{E}_{k}\geq 0$ for any $k$, then $\|\mathcal{E}_{k}\|_{\infty}=\|\mathcal{E}_{k}{\bm{1}}\|_{\infty}$, so that

[TABLE]

Similarly, concerning $F_{2}$ we obtain

[TABLE]

whence

[TABLE]

Concerning $F_{3}$, we have

[TABLE]

whence

[TABLE]

Observe that the reduction of the error per step of the $i$th iteration is bounded from above by $\|H_{i}\|_{\infty}$ which are in turn bounded by $\sigma$ for $i=1,2$ and by $\tau$ for $i=3$. These constants are explicitly computable by means of (15) and (16).

The following result shows that the convergence can be accelerated if $X_{0}$ is a stochastic matrix, say $X_{0}=I$, as in the finite dimensional case [21].

Theorem 4.5**.**

Assume that $A_{-1}{\bm{1}}>A_{1}{\bm{1}}$. For $i\in\{1,2,3\}$ define $X_{k+1}^{(i)}=F_{i}(X_{k}^{(i)})$, $k=0,1,2,\ldots$, where $X_{0}^{(i)}\geq 0$, $X^{(i)}_{0}{\bm{1}}={\bm{1}}$ and $F_{i}(X)$ are given in (14). Then

the three sequences $\{X_{k}^{(i)}\}$ are well defined, 2. 2.

$X_{k}^{(i)}\geq 0$, and $X^{(i)}_{k}{\bm{1}}={\bm{1}}$, 3. 3.

for the error $\mathcal{E}_{k}^{(i)}=G-X_{k}^{(i)}$ we have the following property: $\mathcal{E}_{k}^{(i)}{\bm{1}}=0$, and for any other eigenvector $w\neq{\bm{1}}$ of $G$ such that $Gw=\lambda w$, the sequence $w_{k}^{(i)}=\mathcal{E}_{k}^{(i)}w$ satisfies $\|w_{k+1}^{(i)}\|_{\infty}\leq\|H_{i}(\lambda)\|_{\infty}\|w_{k}^{(i)}\|_{\infty}$, for $i=1,2$ where $H_{1}(\lambda)=(|\lambda|+1)A_{1}+A_{0}$, $H_{2}(\lambda)=(|\lambda|+1)(I-A_{0})^{-1}A_{1}$. Moreover, $\|w_{k+1}^{(3)}\|_{\infty}\leq|\lambda|\cdot\|(I-A_{0}-A_{1}X_{k})^{-1}A_{1}\|_{\infty}\|w_{k}^{(3)}\|_{\infty}$, and $\limsup_{k}\frac{\|w_{k+1}^{(3)}\|_{\infty}}{\|w_{k}^{(3)}\|_{\infty}}\leq|\lambda|$.

Proof 4.6**.**

*We show that if $X\geq 0$ and $X{\bm{1}}={\bm{1}}$ then $F_{i}(X)\geq 0$ and $F_{i}(X){\bm{1}}={\bm{1}}$. For $F_{1}$ this property can be easily checked. For $F_{2}$, since $X\geq 0$ and $(I-A_{0})^{-1}\geq 0$ then $F_{2}(X)\geq 0$. Moreover $F_{2}(X){\bm{1}}=(I-A_{0})^{-1}(A_{-1}+A_{1}){\bm{1}}=(I-A_{0})^{-1}(I-A_{0}){\bm{1}}={\bm{1}}$. Concerning $F_{3}$, since $X{\bm{1}}={\bm{1}}$ then $\|A_{0}+A_{1}X\|_{\infty}=\|A_{0}+A_{1}G\|_{\infty}$ so that in light of Lemma 4.1 the matrix $I-A_{0}-A_{1}X$ is invertible and has nonnegative inverse. This implies that $F_{3}(X)\geq 0$. Moreover, since $X{\bm{1}}={\bm{1}}$ then $(I-A_{0}-A_{1}X-A_{-1}){\bm{1}}=0$ so that $F_{3}(X){\bm{1}}=(I-A_{0}-A_{1}X)^{-1}A_{-1}{\bm{1}}={\bm{1}}$. From (18) we obtain $w_{k+1}^{(1)}=(\lambda A_{1}+A_{1}X_{k}+A_{0})w_{k}^{(1)}$ so that $\|w_{k+1}^{(1)}\|_{\infty}\leq\|\lambda A_{1}+A_{1}X_{k}+A_{0}\|_{\infty}\|w_{k}^{(1)}\|_{\infty}$. On the other hand $\|\lambda A_{1}+A_{1}X_{k}+A_{0}\|_{\infty}\leq\||\lambda|A_{1}+A_{1}X_{k}+A_{0}\|_{\infty}=\|(|\lambda|A_{1}+A_{1}X_{k}+A_{0}){\bm{1}}\|_{\infty}=\|(|\lambda|A_{1}+A_{1}+A_{0}){\bm{1}}\|_{\infty}=\|H_{1}(\lambda)\|_{\infty}$. Similarly, we proceed with $F_{2}$ relying on (19). Concerning $F_{3}$, from (20) we have $w_{k+1}^{(3)}=\lambda(I-A_{0}-A_{1}X_{k})^{-1}A_{1}w_{k}^{(3)}$, whence $\|w_{k+1}^{(3)}\|_{\infty}\leq|\lambda|\|(I-A_{0}-A_{1}X_{k})^{-1}A_{1}\|_{\infty}\|w_{k}^{(3)}\|_{\infty}$. Since $A_{1}{\bm{1}}<A_{-1}{\bm{1}}$ then $\|(I-A_{0}-A_{1}X_{k})^{-1}A_{1}\|_{\infty}=\|(I-A_{0}-A_{1}X_{k})^{-1}A_{1}{\bm{1}}\|_{\infty}\leq\|(I-A_{0}-A_{1}X_{k})^{-1}A_{-1}{\bm{1}}\|_{\infty}=1$. Taking the limsup for $k\to\infty$ we obtain $\limsup_{k}\|(I-A_{0}-A_{1}X_{k})^{-1}A_{1}\|_{\infty}\leq\|(I-A_{0}-A_{1}G)^{-1}A_{-1}\|_{\infty}=\|G\|_{\infty}=1$. This completes the proof. *

Observe that the condition $\lim_{k}\|G-X_{k}\|_{\infty}=0$ implies that $\|w_{k+1}^{(3)}\|_{\infty}/\|w_{k}^{(3)}\|_{\infty}\leq|\lambda|\lim_{k}\|(I-A_{0}-A_{1}X_{k})^{-1}A_{1}\|_{\infty}=|\lambda|\cdot\|H_{3}\|_{\infty}\leq|\lambda|\cdot\tau$.

According to the above theorem, the sequences generated by the three functional iterations with $X_{0}$ stochastic, converge faster than the corresponding sequences obtained with $X_{0}=0$. For the functions $F_{1}$ and $F_{2}$, this follows since $\|H_{1}(\lambda)\|_{\infty}\leq\|H_{1}\|_{\infty}$ and $\|H_{2}(\lambda)\|_{\infty}\leq\|H_{2}\|_{\infty}$. Similarly, we may proceed for the function $F_{3}$. The convergence of the sequence is faster the smaller $\sup\{|\lambda|:\,\lambda\neq 1,~{}\lambda\hbox{ ei}\hbox{gen}\hbox{va}\hbox{lue of }G\}$.

4.1 Implementation issues

Let $g(z)$ be the solution of minimum modulus of (7) and consider the sequences generated by the fixed point iterations (14) obtained starting with $X_{0}=T(g)+C$, where $C$ is any correction. Denote by $X_{k}$ any one of these three sequences so that we have $X_{k}=T(g)+E_{k}$, $E_{0}=C$, where $E_{k}$ is the correction part. In this section we aim to explicit the equation which relates $E_{k+1}$ to $E_{k}$.

Consider the first iteration $X_{k+1}=A_{1}X_{k}^{2}+A_{0}X_{k}+A_{-1}$ and denote by $F$ the correction matrix such that $T(g)+F=A_{1}T(g)^{2}+A_{0}T(g)+A_{-1}$. Subtracting the latter equation from the former and performing formal manipulations yields

[TABLE]

This equation provides a more efficient way to implement the first iteration since it involves multiplications of QT matrices and correction matrices and avoids the multiplication of QT matrices having a nonzero symbol. In this version, the precomputation of $g$, $S$ and $F$ is needed.

Consider the second iteration $X_{k+1}=(I-A_{0})^{-1}(A_{1}X_{k}^{2}+A_{-1})$. Subtract it from the equation $T(g)=(I-A_{0})^{-1}(A_{1}T(g)^{2}+A_{-1})-(I-A_{0})^{-1}F$ and performing formal manipulations yields

[TABLE]

Also in this case the precomputation of $T(g)$, $\widehat{S}$, $F$ and $\widetilde{S}$ is needed.

Concerning the third iteration and proceeding similarly we arrive at the recursion

[TABLE]

Also in this case the precomputation of $T(g)$, $\widehat{V}$, and $\widetilde{V}$ is needed. However, at each step the inverse of a QT matrix must be computed.

The iterations (21)–(23) can be started with $X_{0}=T(g)$, that is, $E_{0}=0$. Alternatively, they can be started with $X_{0}=T(g)+ve_{1}^{T}$, where $e_{1}^{T}=(1,0,\ldots)$ and $v$ is chosen in such a way that $X_{0}$ is row-stochastic so that convergence is faster. This can be accomplished by setting $E_{0}=ve_{1}^{T}$.

5 Newton’s method

Rewrite equation $A_{1}X^{2}+A_{0}X+A_{-1}=X$ as

[TABLE]

Newton’s method applied to equation (24) generates the sequence

[TABLE]

where the matrix $Z_{k}$ solves the equation $L^{\prime}(Z_{k})=L(X_{k})$, and $L^{\prime}(H)=A_{1}XH+A_{1}HX+(A_{0}-I)H$ is the Fréchet derivative of $L(X)$ at $X$ applied to the matrix $H$. More specifically, $Z_{k}$ solves the following Sylvester equation

[TABLE]

Observe that we may write

[TABLE]

provided that $\|(I-A_{0}-A_{1}X_{k})^{-1}\|_{\infty}$ is bounded from above, $\|S_{k}\|_{\infty}<1$ and $\|X_{k}\|_{\infty}\leq 1$ so that we have

[TABLE]

Moreover, we have

[TABLE]

The latter equality can be proved by observing that in general, for $Y=X-H$, we have $L(Y)=A_{1}(X-H)^{2}+(A_{0}-I)(X-H)+A_{-1}=L(X)-L^{\prime}(H)+A_{1}H^{2}$ so that, if $H$ is such that $L^{\prime}(H)=L(X)$ as in a Newton step, then $L(Y)=A_{1}H^{2}$.

Another useful property is the following. Equation (25) can be rewritten as $Z_{k}=\mathcal{E}_{k+1}-\mathcal{E}_{k}$, where $\mathcal{E}_{k}=G-X_{k}$, and $L(X_{k})$ can be rewritten as $L(X_{k})=L(X_{k})-L(G)=-A_{1}(\mathcal{E}_{k}G+X_{k}\mathcal{E}_{k})-(A_{0}-I)\mathcal{E}_{k}$. Replace these two representations for $Z_{k}$ and $L(X_{k})$ in (26), and get

[TABLE]

By following the same arguments used to arrive at (27), we may rewrite the above equation as

[TABLE]

and get

[TABLE]

The following result extends to QT matrix coefficients the convergence results valid in the finite dimensional case [16]:

Theorem 5.1**.**

Assume that $A_{-1}{\bm{1}}>A_{1}{\bm{1}}$. Let $X_{k}$, $k=0,1,\ldots,$ be the sequence generated by (25) and (26) starting with $X_{0}=0$. Then, for any $k=0,1,2,\ldots$,

equation (26) has a solution $Z_{k}$ such that $\|Z_{k}\|_{\infty}\leq\beta$, where

[TABLE]

*for $W=I-A_{0}-A_{1}G$, so that $X_{k+1}$ is well defined; * 2. 2.

$Z_{k}\leq 0$, $L(X_{k+1})\geq 0$ and $0\leq X_{k}\leq X_{k+1}\leq G$; 3. 3.

$\lim_{k\to\infty}(\mathcal{E}_{k})_{i,j}=0$* for any $i,j\geq 1$, where $\mathcal{E}_{k}=G-X_{k}$;* 4. 4.

$\|\mathcal{E}_{k+1}\|_{\infty}\leq\frac{\tau}{1-\tau}\|\mathcal{E}_{k}\|^{2}_{\infty}$, $\|Z_{k+1}\|_{\infty}\leq\frac{\tau}{1-\tau}\|Z_{k}\|^{2}_{\infty}$, where $\tau=\min\{\gamma,\sigma\}$, with $\gamma$ and $\sigma$ defined in (6) and (15), respectively.

Proof 5.2**.**

We prove properties 1 and 2 by induction on $k$. For $k=0$ we have $X_{0}=0$, $Z_{0}=-(I-A_{0})^{-1}A_{-1}\leq 0$, $L(X_{1})=A_{1}Z_{0}^{2}\geq 0$, $X_{1}=-Z_{0}\geq X_{0}$ and $X_{1}=(I-A_{0})^{-1}A_{-1}\leq G$. Moreover, clearly $\|Z_{0}\|_{\infty}\leq\beta$. For the inductive step, assume that properties 1 and 2 are valid for $k$ and prove them for $k+1$. We show that $\|Z_{k+1}\|\leq\beta$. Consider (28). Since by induction $L(X_{k+1})\geq 0$ then $\|L(X_{k+1})\|_{\infty}=\|L(X_{k+1}){\bm{1}}\|_{\infty}\leq\|(A_{1}X_{k+1}^{2}+A_{0}X_{k+1}+A_{-1}){\bm{1}}\|_{\infty}+\|X_{k+1}{\bm{1}}\|_{\infty}$. Moreover, since $X_{k+1}\leq G$ then $X_{k+1}{\bm{1}}\leq G{\bm{1}}$ so that $\|L(X_{k+1})\|_{\infty}\leq 2$. From $X_{k+1}\leq G$ it follows also $A_{0}+A_{1}X_{k+1}\leq A_{0}+A_{1}G$ so that $(I-A_{0}-A_{1}X_{k+1})^{-1}\leq(I-A_{0}-A_{1}G)^{-1}=W^{-1}$, whence $\|(I-A_{0}-A_{1}X_{k+1})^{-1}\|_{\infty}\leq\|W^{-1}\|_{\infty}$. Similarly, $\|S_{k}\|_{\infty}\leq\|(I-A_{0}-A_{1}G)^{-1}A_{1}\|_{\infty}=\|W^{-1}A_{1}\|_{\infty}<1$ in view of Lemma 4.1. From (27) and (28) we get $\|Z_{k}\|_{\infty}\leq\beta$. The property $Z_{k+1}\leq 0$ follows from (27) since $S_{k+1}\geq 0$, $(I-A_{0}-A_{1}X_{k+1})^{-1}\geq 0$, and $L(X_{k+1})\geq 0$ by the inductive assumption. The inequality $L(X_{k+2})\geq 0$ follows from (29) since $Z_{k+1}\leq 0$. Consequently $X_{k+2}=X_{k+1}-Z_{k+1}\geq X_{k+1}$. The property $X_{k+2}\leq G$ follows from (30) since $\mathcal{E}_{k+1}\geq 0$, $S_{k+1}\geq 0$ and $X_{k+1}\geq 0$. Given $i,j$ consider the sequence $(\mathcal{E}_{k})_{i,j}$ for $k=0,1,\ldots$. This sequence is non-increasing and bounded from below by [math] therefore it has a limit. The value of the limit cannot be positive since $G$ is the minimal nonnegative solution to the matrix equation. From the representation (30) of $\mathcal{E}_{k+1}$, since $\|X_{k}\|_{\infty}\leq\|G\|_{\infty}\leq 1$, and $\|S_{k}\|_{\infty}\leq\|(I-A_{0}-A_{1}G)^{-1}A_{1}\|_{\infty}\leq\tau<1$ in view of (16), we deduce that

[TABLE]

*Similarly we can do for $\|Z_{k+1}\|_{\infty}$. *

6 Perturbation results

For the case where the coefficient matrices are finite, Higham and Kim [13] derived a condition number $\Psi(X)$ for a solvent $X$ of a general quadratic matrix equation of the kind (1) namely,

[TABLE]

where $P=I_{n}\otimes A_{1}X+X^{T}\otimes A_{1}+I_{n}\otimes(A_{0}-I)$ and $\alpha,\beta,\gamma$ are nonnegative parameters.

However, when the coefficient matrices are semi-infinite, there are cases where $P^{-1}$ does not exist or $\|X\|_{F}=\infty$, so the definition of $\Psi(X)$ does not apply. In this section, we take into account the structure of the coefficient matrices and derive a structured condition number for the minimal nonnegative solution of equations (1) and (2). Without loss of generality we consider only equation (1).

Consider the perturbed matrix equation obtained from (1) by replacing the coefficients $A_{i}$ by $A_{i}+\Delta_{A_{i}}$ where $\Delta_{A_{i}}=T(\delta_{i})+E_{\delta_{i}}\in\mathcal{QT},$ $A_{i}+\Delta_{A_{i}}\geq 0$, for $i=-1,0,1$, and $(A_{1}+\Delta_{A_{1}}+A_{0}+\Delta_{A_{0}}+A_{-1}+\Delta_{A_{-1}})\bf{1}=\bf{1}$. Denote $X+\Delta_{X}$ a solution of the perturbed equation so that we may write

[TABLE]

The analysis is separated into two parts, that is, the the analysis of the structured condition number of the Toeplitz part and the analysis of the condition number of the whole matrix.

6.1 Toeplitz part

In this section we provide a perturbation result for the function $g(z)$ which is the solution of minimum modulus of the scalar equation (7).

For the sake of notational simplicity, we omit the variable $z$ from the symbols, say, we write $g$ in place of $g(z)$ and $a_{i}$ in place of $a_{i}(z)$.

Under the assumption that the matrix coefficients $A_{i}+\Delta_{A_{i}}$ of equation (31) still satisfy the condition $A_{i}+\Delta_{A_{i}}\geq 0$, $\sum_{i=-1}^{1}(A_{i}+\Delta_{A_{i}}){\bm{1}}={\bm{1}}$, for Theorem 2.3 the minimal nonnegative solution of (31) can be written as $G+\Delta_{G}$ , where $\Delta_{G}=T(\delta_{g})+E_{\delta_{g}}\in\mathcal{QT}$ and $g+\delta_{g}$ is the solution of minimum modulus of the equation

[TABLE]

Taking the difference of the above equation with (7), where we set $\mu=g+\delta_{g}$ and $\lambda=g$, we obtain

[TABLE]

Whence, neglecting higher order terms in the perturbations we get

[TABLE]

where $\doteq$ means equality up to higher order terms with respect to the perturbations. This yields

[TABLE]

Note that $g(z)$, $a_{0}(z)$ and $a_{1}(z)$ have nonnegative coefficients so that $\|g\|_{w}=\sum_{i\in\mathbb{Z}}g_{i}=g(1)=1$, and $|2a_{1}(z)g(z)+a_{0}(z)|\leq 2a_{1}(1)g(1)+a_{0}(1)=2a_{1}(1)+a_{0}(1)=1-a_{-1}(1)+a_{1}(1)<1$, due to (8). This way, we have $\|(1-2a_{1}g-a_{0})^{-1}\|_{w}=(1-2a_{1}(1)g(1)-a_{0}(1))^{-1}=(a_{-1}(1)-a_{1}(1))^{-1}$. Whence, from (32), we obtain

[TABLE]

where $\operatorname*{\,\dot{\leq}\,}$ means inequality up to higher order terms with respect to the perturbations. Therefore we arrive at the bound

[TABLE]

If we measure the perturbations by

[TABLE]

we have $\|\delta_{i}\|_{w}\leq\epsilon\|a_{i}\|_{w}$, moreover, since $\|a_{-1}\|_{w}+\|a_{0}\|_{w}+\|a_{1}\|_{w}=a_{-1}(1)+a_{0}(1)+a_{1}(1)=1$, the relative variation of the symbol is bounded by

[TABLE]

It follows from (34) that $\hbox{cond}_{T(g)}:=\frac{1}{a_{-1}(1)-a_{1}(1)}$ is an upper bound to the condition number of the Toeplitz part of $G$.

6.2 Whole matrix

Expanding (31), omitting the second and higher order terms in the perturbations, and setting $X=G$ lead to

[TABLE]

Now we prove some properties that will be useful to estimate the condition number of the whole solution.

Let $\Delta_{A}:=\Delta_{A_{1}}G^{2}+\Delta_{A_{0}}G+\Delta_{A_{-1}}$. According to (5), equation (35) can be written as

[TABLE]

where $F:\mathcal{QT}\rightarrow\mathcal{QT}$ is the map defined by

[TABLE]

Now, we prove that the map $F(Y)$ is invertible in $\mathcal{L}_{\infty}$, that is, $F^{-1}$ has bounded infinity norm. By Lemma 2.1, we have $\|W^{-1}RW\|_{\infty}\leq\gamma<1$ and $\|G\|_{\infty}=1$ so that the series $\sum_{k=0}^{\infty}(W^{-1}RW)^{k}VG^{k}$ is convergent for any $V\in\mathcal{L}_{\infty}$. Therefore, if $V=F(Y)$ then $Y=\sum_{k=0}^{\infty}(W^{-1}RW)^{k}VG^{k}$. Thus, we get

[TABLE]

Since $\sum_{k=0}^{\infty}\|W^{-1}RW\|_{\infty}^{k}=1/(1-\|W^{-1}RW\|_{\infty})$, taking norms in the above expression and applying Lemma 2.1 yields

[TABLE]

Whence we conclude with the following

Theorem 6.1**.**

If $A_{-1}{\bm{1}}>A_{1}{\bm{1}}$, then for the perturbation $\Delta_{G}$ we have

[TABLE]

*where $\theta$ and $\gamma$ are defined in (6) and $\Delta_{A}=\Delta_{A_{1}}G^{2}+\Delta_{A_{0}}G+\Delta_{A_{-1}}$. *

From the above result it turns out that $\|W^{-1}\|_{\infty}/(1-\|W^{-1}RW\|_{\infty})$ is an estimate of the conditioning of the problem, while $1/(\theta(1-\gamma))$ provides an upper bound. Since $1-a_{1}(1)/a_{-1}(1)=(a_{-1}(1)-a_{1}(1))/a_{-1}(1)$, we may rewrite the upper bound to the conditioning in the following form which is closer to the expression obtained for the Toeplitz part of $G$ in Section 6.1.

[TABLE]

It is interesting to observe that if $a_{-1}(1)\leq b_{-1}(1)$ and $\frac{a_{-1}(1)}{b_{-1}(1)}\leq\frac{a_{1}(1)}{b_{1}(1)}$, which in turn is verified if $a_{-1}(1)\leq b_{-1}(1)$ and $b_{1}(1)\leq a_{1}(1)$, then

[TABLE]

that is, the conditioning of the Toeplitz part coincides with the conditioning of the whole problem.

It is interesting to point out that, since $RW=A_{1}$ (see equation (5)), the estimate of the condition number $\|W^{-1}\|_{\infty}/(1-\|W^{-1}RW\|_{\infty})$, appears in the uniform bound $\beta$ to the norm of the Newton correction $Z_{k}$ introduced in Theorem 5.1. Consequently, if the quadratic matrix equation is well conditioned, the uniform bound to the norm of $Z_{k}$ is smaller.

6.3 A simple example

This example is taken from Example 6.2 in [24], where a continuous time Markov process modeling a two-node Jackson network is considered. Here, the matrices are modified by means of the uniformization technique [19] in order to represent a discrete-time model. In details,

[TABLE]

[TABLE]

where the parameters $\lambda_{1},\lambda_{2},\mu_{1},\mu_{2},p,q$ are chosen as in Table 1 and $\alpha=(\lambda_{1}+\lambda_{2}+\mu_{1}+\mu_{2})^{-1}$. Each case denotes one instance of the two-node Jackson network, formed by two servers and two queues, where customers arrive at nodes 1 and 2 according to independent Poisson processes with rates $\lambda_{1}$ and $\lambda_{2}$, respectively. Customers are served according to a first-come-first-served discipline, service times at nodes 1 and 2 are independent and exponentially distributed with means $1/\mu_{1}$, $1/\mu_{2}$. After completing service at node 1, customers enter node 2 with probability $p$ or leave the system with probability $1-p$, where $0<p<1$. After completing service at node 2, customers enter node 1 with probability $q$ or leave the system with probability $1-q$, where $0<q<1$.

In [24], 10 cases defined by the parameters given in Table 1 are analyzed. It can be easily seen that the condition $A_{-1}{\bm{1}}>A_{1}{\bm{1}}$ holds for cases 1, 3, 4, 5, 7, 8, 9, while the same condition holds in the cases 2, 6 and 10 for the flipped problems where phases and levels are exchanged.

For the flipped problem, the coefficient matrices are

[TABLE]

[TABLE]

If we perturb the parameters $\lambda_{1},\lambda_{2},\mu_{1},\mu_{2}$ by

[TABLE]

where the perturbations are randomly chosen, then we get the perturbed matrices $A_{-1}+\Delta_{A_{-1}}$, $A_{0}+\Delta_{A_{0}}$ and $A_{1}+\Delta_{A_{1}}$.

Note that $a_{-1}(1)=b_{-1}(1)$, $a_{-1}(1)-a_{1}(1)<b_{-1}(1)-b_{1}(1)$ holds true for both the original and the flipped problems, it follows $\frac{1}{\theta(1-\gamma)}=\frac{1}{a_{-1}(1)-a_{1}(1)}$, that is, the conditioning of the Toeplitz part coincides with the conditioning of the whole problem.

We denote by “Conditioning” the upper bound on the condition number for the minimal nonnegative solution $G$, that is, $\frac{1}{\theta(1-\gamma)}$. Moreover, we denote by $\delta_{g}$-bound and $\Delta_{G}$-bound, respectively, the perturbation bound (33) on $\|\delta_{g}\|_{w}$ and the bound (38) on $\|\Delta_{G}\|_{\infty}$. In Table 2, we report the upper bound on the condition number for the minimal nonnegative solution $G$ of equation (1), and we compare the perturbation bound (33) on the Toeplitz part of $G$ and the bound on the solution $G$ with the corresponding perturbation errors.

It can be seen from Table 2 that the upper bound $\frac{1}{\theta(1-\gamma)}$ can serve as a very good estimate of the conditioning of the problem. Inequalities (33) and (38) provide very sharp and revealing perturbation bounds to the Toeplitz part and to the solution $G$ with respect to small perturbations on the coefficients.

7 Computational issues and numerical experiments

Observe that, since the class $\mathcal{QT}$ is an algebra, then all the matrices $X_{k}^{(i)}$ generated by the fixed point iterations of Section 4 belong to $\mathcal{QT}$ and each fixed point iteration can be easily implemented in Matlab relying on the CQT-Toolbox of [6]. Concerning Newton iteration, a crucial role is played by the solution of the Sylvester-like equation (26). In the case where the coefficients have finite size $n$, this equation can be solved by the Bartels-Stewart algorithm [1] in $O(n^{3})$ arithmetic operations. The case of infinite size coefficients is much more complicated. For quasi-Toeplitz matrices, the problem has been analyzed in [28] by using rational Krylov subspaces techniques where it is proved that $Z_{k}\in\mathcal{QT}$ under suitable assumptions that are satisfied by the condition $A_{-1}{\bm{1}}>A_{1}{\bm{1}}$. In the numerical experiments reported in this section we have used the implementation of [28] to solve (26).

We have implemented the three fixed point iterations analyzed in Section 4 both in the standard versions (14), and in the version which separately computes the correction part, see (21), (22) and (23). We also implemented the computation of the coefficients of $g(z)$ relying on the evaluation and interpolation strategy at the roots of 1 with automatic handling of the number of interpolation points described in Algorithm 1.

For each fixed point iteration, we have considered four different starting approximations, namely, $X_{0}=0$, $X_{0}=I$, $X_{0}=T(g)$, $X_{0}=T(g)+ve_{1}^{T}$ where $e_{1}=(1,0,\ldots)^{T}$ and $v$ is chosen so that $X_{0}$ is row-stochastic.

Since there is not much difference in the performances of the versions based on computing only the correction and the versions where the whole matrix is computed, we report only the results concerning the latter versions.

We have compared the three fixed point iterations, with different choices of $X_{0}$, to Newton’s iteration and to the algorithm of cyclic reduction (CR) analyzed in [4], which, in the case of finite matrices, is the method of choice commonly used in practice. For each experiment, we report the number of iterations, and CPU time needed to reach the bound $\|A_{1}X^{2}+(A_{0}-I)X+A_{-1}\|_{\infty}\leq\epsilon$ to the residual error where $\epsilon=5.0\cdot 10^{-14}$. We have considered some test problems modeling real world networks. More precisely, the “Two-node Jackson network” of Example 6.2 in [24], reported in Section 6.3, and the model “Assistance from an idle server” of [24] with different choices of the parameters, together with a general random walk in the quarter plane where the assigned probabilities have been chosen in such a way to have long queues in the system.

7.1 Two-node Jackson network

The model that we consider has been described in Section 6.3. Among the 10 problems in the list of Table 1, we report the results of Problem 7 which is the most ill-conditioned in the list, together with the case obtained with different values of the parameters $\lambda_{i},\mu_{i}$, $p$ and $q$ which make the matrix $G$ numerically very large so that the computational effort is substantially large.

Figure 1 concerns Problem 7 in the list of Table 1. The first three graphs report the residual errors at each step for different values of $X_{0}$. The fourth graph compares the residual errors of the three iterations for $X_{0}=T(g)+ve_{1}^{T}$. In Table 3 it is reported the CPU time in seconds (left) together with the number of steps (right) required by the three iterations to arrive at a residual error at most $5.0\cdot 10^{-14}$. The computation of the symbol $g(z)$ is very inexpensive since the coefficients $g_{i}$ are computed in $0.003$ seconds. For this problem, cyclic reduction provides the solution in just 8 steps and in 1.97 seconds, while Newton iteration requires 8 steps but takes a larger amount of seconds, i.e., 256.8.

We may observe that in this case CR is the most efficient algorithm and that the results of Theorems 4.3 and 4.5 are respected. In fact the iteration given by $F_{3}$ with $X_{0}=T(g)+ve_{1}^{T}$ is the fastest one in terms of number of steps. While concerning the CPU time, the iteration $F_{2}$ with $X_{0}=T(g)+ve_{1}^{T}$ is the fastest.

Table 4 concerns the case of a Two-node Jackson network with the choice of parameters given by $\lambda_{1}=5$, $\lambda_{2}=0.7$, $\mu_{1}=2$, $\mu_{2}=2$, $p=0.5$ and $q=0.5$. In this model, seen as a random walk in the quarter plane, the overall probability to move right is higher than the overall probability to move left, while the probability to move down is higher than the probability to move up. The table reports the CPU time in seconds (left) together with the number of steps required by the three iterations (right) to arrive at a residual error at most $5.0\cdot 10^{-14}$. The computation of the symbol $g(z)$ remains almost inexpensive even though the numerical length of the coefficient vector of the symbol $g(z)$ is quite large, in fact the coefficients $g_{i}$ are computed in less than $0.007$ seconds, and the size of the coefficient vector is $31$ for the coefficients of the negative powers of $z$ and $8424$ for the coefficients of the positive powers, respectively. The numerical size of the correction is $28\times 6937$.

For this problem, cyclic reduction provides the solution in 8 steps and in $70.8$ seconds, while Newton iteration requires 8 steps and takes 782 seconds. In this case, due to the large size of the matrices involved in the computation of matrix inverses, CR takes a much larger time than the simple functional iterations given by $F_{1}$ and $F_{2}$ which either do not involve inversion or require just only one matrix inversion. While iteration given by $F_{2}$ takes less than 10 seconds, the iteration given by $F_{3}$, even though is the fastest in terms of number of steps, needs a large CPU time. In fact, similarly to CR, it requires a matrix inversion at each step which becomes more expensive as the approximation approaches the limit. Newton iteration has the same convergence features as CR, however, the larger cost of solving a Sylvester equation makes this iteration much slower than the other ones at least for this problem.

In this model, increasing the values of $\lambda_{1}$ makes the size of the output much larger. In particular, with values $\lambda_{1}\geq 6$, cyclic reduction breaks down for memory overflow while $F_{2}(X)$ provides the solution with a slight increase of the CPU time.

Figure 2 provides the solution $G$ in log scale where the Toeplitz part and the correction parts are separately represented.

7.2 Assistance from idle server

Here we consider a class of queueing models for a system with two servers and two queues. Arrivals to queues 1 and 2 occur as independent Poisson processes with parameters $\lambda_{1}$ and $\lambda_{2}$, respectively. The service times of servers 1 and 2 are exponentially distributed with parameters $\mu_{1}$ and $\mu_{2}$ respectively. Each server serves its own queue according to a first-come-first-served discipline. If one of the queues is empty, the server for that queue assists the other server, doubling the latter’s service rate. If there is an arrival to a queue while its server is assisting the other queue, the server immediately ceases assisting and serves its own queue. This stochastic process is ergodic if and only if $\rho_{1}+\rho_{2}<2$, where $\rho_{i}=\lambda_{i}/\mu_{i}$, $i=1,2$.

For this model, the matrices $A_{-1},A_{0},A_{1}$ are given by $A_{-1}=\hbox{diag}(2\mu_{1},\mu_{1},\mu_{1},\ldots)$, $A_{0}=\hbox{trid}(\mu_{2},-\lambda_{1}-\lambda_{2}-\mu_{1}-\mu_{2},\lambda_{2})+(\mu_{2}-\mu_{1})e_{1}e_{1}^{T}$, $A_{1}=\lambda_{1}I$.

In this example, we have chosen the values of the parameters in order that the numerical size of the matrix $G$ is substantially large, namely, $\lambda_{1}=0.01$, $\lambda_{2}=2.9$, $\mu_{1}=0.03$, $\mu_{2}=2.0$.

The situation is analogous to that of the second example in Section 7.1 where the methods based on functional iterations perform better than cyclic reduction and Newton iteration. The graphs in Figure 3 and the values in Table 5 synthesize the behavior of the algorithms. Cyclic reduction takes about 17 seconds of CPU while Newton iteration about 600 seconds. Slightly increasing the value of $\lambda_{1}$, CR breaks down for memory overflow while functional iterations still compute correctly the solution $G$.

7.3 Random walk in the quarter plane

Here we consider an example where the condition $A_{-1}{\bm{1}}>A_{1}{\bm{1}}$ is satisfied everywhere exept in the first component. The example describes a random walk in the quarter plane where a particle can occupy positions in a grid and we know the probabilities that the particle moves to the neighboring positions. In order to better describe the test problem, we denote $H=(h_{i,j})_{i,j=-1,1}$ the matrix with the probabilities of transition in the inner part of the quarter plane, while we denote $Y=(y_{i,j})$ the $3\times 2$ matrix with the probabilities of transition in the $y$ axis. These two matrices fully describe the coefficients $A_{i}$ which can be written as $A_{i}=T(a_{i})+E_{i}$ where $a_{i}(z)=\sum_{j=-1}^{1}h_{i,j}z^{j}$ and $E_{i}=e_{1}[y_{i,0}-h_{i,0},y_{i,1}-h_{i,1},0,\ldots]$, compare with (4).

The random walk of this example is obtained with the values

[TABLE]

so that the condition $A_{-1}{\bm{1}}>A_{1}{\bm{1}}$ is satisfied in all the components but the first.

For this problem, there exist two nonnegative solutions $G$ and $\widehat{G}$ to equation (1) which satisfy the inequality $G\leq\widehat{G}$. Moreover $\widehat{G}$ is stochastic while $G$ is substochastic. The two solutions have the same symbol $g(z)$ and differ only for the correction part. Starting with $X_{0}=0$ or with $X_{0}=T(g)$ the sequences generated by the functional iterations converge to $G$. Starting with $X_{0}=I$ or $X_{0}=T(g)+ve_{1}^{T}$ the sequences converge to $\widehat{G}$, while CR and Newton iteration converge to $G$.

The results of this test are summarized in Table 6. Both CR and Newton iteration take 23 steps in order to arrive at numerical convergence. However CR takes more than 8 minutes of CPU time while Newton iteration just 3.4 seconds. The large amount of CPU time taken by CR is due to the fact that the inverse matrices involved at each step of CR are QT matrices with a correction having a size which increases step after step and reaches values larger than $10^{6}$, while Newton iteration involves QT matrices with corrections having almost the same size of the correction of $G$ which is $126\times 36$. The growth of the sizes of the correction matrices in the algorithms is an issue which deserves further analysis.

Functional iterations $F_{1},F_{2}$ and $F_{3}$ with $X_{0}=0$ or with $X_{0}=T(g)$ take a large number of steps to converge numerically to $G$ while with $X_{0}=I$ or $X_{0}=T(g)+ve_{1}^{T}$ the number of iterations is much smaller but the limit of the sequences is the stochastic solution $\widehat{G}$ which is not the minimal one.

For this problem, the combination of few steps of CR followed by few steps of Newton iteration provide a substantial acceleration in terms of CPU time. In fact, the first iterations of CR, involving matrices of small size, have a low cost. The last few steps of CR, which have a much higher cost, are replaced by Newton steps.

8 Conclusions

We have analyzed quadratic matrix equations encountered in the solution of random walk in the quarter plane where the solution of interest is the minimal nonnegative solution $G$. This class of equations is characterized by matrix coefficients with infinite size which belong to the class $\mathcal{QT}$ of Quasi-Toeplitz matrices. We have provided a perturbation analysis of $G$, introduced some fixed point algorithms for computing $G$ and compared their convergence speed. The algorithms rely on the properties of $\mathcal{QT}$ matrices recently investigated in [6]. Numerical experiments show that in many cases the CPU time and the memory resources required by our approach are significantly inferior to the ones required by the algorithm of cyclic reduction, which is considered as the algorithm of choice for this class of problems. The effectiveness of Newton iteration depends on the growth of the sizes of the correction part in the QT matrices generated by the algorithm.