Relaxed Fixed Point Iterations for Matrix Equations Arising in Markov Chains Modeling^††thanks: The authors are partially supported by INDAM/GNCS and by the project PRA_2020_61 of the University of Pisa.

The transition probability matrix of an M/G/1-type Markov chain is a block Hessenberg matrix $P$ of the form

with $A_{i},B_{i}\in\mathbb{R}^{n\times n}\geq 0$ and $\sum_{i=0}^{\infty}B_{i}$ and $\sum_{i=-1}^{\infty}A_{i}$ stochastic matrices.

In the sequel, given a real matrix $A=(a_{ij})_{ij}\in\mathbb{R}^{m\times n}$, we write $A\geq 0$ ($A>0$) if $a_{ij}\geq 0$ ($a_{i,j}>0$) for any $i,j$. A stochastic matrix is matrix $A\geq 0$ such that $A\mbox{\boldmath$e$}=\mbox{\boldmath$e$}$, where $e$ is the column vector having all the entries equal to 1.

In the positive recurrent case, the computation of the steady state vector $\pi$ of $P$, such that

is related with the solution of the unilateral power series matrix equation

Under some mild assumptions this equation has a componentwise minimal non-negative solution $G$ which determines, by means of Ramaswami’s formula [16], the vector $\pi$.

Among the easy-to-use, but still effective, tools for numerically solving (3), there are fixed point iterations (see [3] and the references given therein for a general review of these methods). The intrinsic simplicity of such schemes make them attractive in domains where high performance computing is crucial. But they come at a price: the convergence can become very slow especially for problems which are close to singularity. The design of acceleration methods (aka extrapolation methods) for fixed point iterations is a classical topic in numerical analysis [5]. Relaxation techniques are commonly used for the acceleration of classical stationary iterative solvers for large linear systems. In this paper we introduce some new coupled fixed point iterations for solving (3) which can be combined with relaxation techniques to speed up their convergence.

More specifically, we first observe that computing the solution of the matrix equation (3) is formally equivalent to solving a semi-infinite block Toeplitz block Hessenberg linear system. Customary block iterative algorithms applied for the solution of this system yield classical fixed point iterations. In particular, the traditional and the U-based fixed point iteration [3] originate from the block Jacobi and the block Gauss-Seidel method, respectively. Recently, in [7] the authors show that some iterative solvers based on a block staircase partitioning outperform the block Jacobi and the block Gauss-Seidel method for M-matrix linear systems in block Hessenberg form. The application of the staircase splitting to the block Toeplitz block Hessenberg linear system associated with (3) yields the new coupled fixed point iteration

starting from an initial approximation $X_{0}$. The contribution of this paper is aimed at highlighting the properties of (4).

We show that, if $X_{0}=0$, the sequence $\{X_{k}\}_{k}$ defined in (4) converges to $G$ faster than the traditional fixed point iteration. In the case where the starting matrix $X_{0}$ of (4) is any column stochastic matrix and $G$ is also column stochastic, we prove that the sequence $\{X_{k}\}_{k}$ still converges to $G$. Moreover, by comparison of the mean asymptotic rates of convergence, we conclude that (4) is asymptotically faster than the traditional fixed point iteration.

At each iteration the scheme (4) determines two approximations which can be combined by using the relaxation technique, that is, the approximation computed at the $k$-th step takes the form of a weighted average between $Y_{k}$ and $X_{k+1}$. The modified relaxed variant is defined by the sequence

where $\omega_{k+1}$ is the relaxation parameter. The convergence results for (4) easily extend to the modified scheme in the case of under-relaxation, that is, the parameter $\omega_{k}$ satisfies $0\leq\omega_{k}\leq 1$. Heuristically, it is argued that over-relaxation values ($\omega_{k}>1)$ can improve the convergence. If $X_{0}=0$, under some assumptions a theoretical estimate of the asymptotic convergence rate of (5) is given which confirms this heuristic. Moreover, an adaptive strategy is devised which makes possible to perform over-relaxed iterations of (5) by still ensuring the convergence of the overall iterative process. The results of extensive numerical experiments confirm the effectiveness of the proposed variants, which generally outperform the $U$-based fixed point iteration for nearly singular problems. In particular, the over-relaxed scheme (5) with $X_{0}=0$, combined with the adaptive strategy for parameter estimation, is capable to significantly accelerate the convergence without increasing the computational cost.

The paper is organized as follows. In Section 2 we set up the theoretical framework, briefly recalling some preliminary properties and assumptions. In Section 3 we revisit classical fixed point iterations for solving (3) by establishing the link with the iterative solution of an associated block Toeplitz block Hessenberg linear system. In Section 4 we introduce the new fixed point iteration (4) by proving some convergence results. The relaxed variant (5), as well as the generalizations of convergence results for this variant, are described in Section 5. Section 6 deals with a formal analysis of the asymptotic convergence rate of both (4) and (5). Adaptive strategies for the choice of the relaxation parameter are discussed in Section 7, together with their cost analysis under some simplified assumptions. Finally, the results of extensive numerical experiments are presented in Section 8 whereas conclusions and future work are the subjects of Section 9.

Throughout the paper we assume that $A_{i}$, $i\geq-1$, are $n\times n$ nonnegative matrices, such that their sum $A=\sum_{i=-1}^{\infty}A_{i}$ is irreducible and row stochastic, that is, $A\mbox{\boldmath$e$}=\mbox{\boldmath$e$}$, $\mbox{\boldmath$e$}=\left[1,\ldots,1\right]^{T}$. According to the results of [3, Chapter 4], such assumption implies that (3) has a unique componentwise minimal nonnegative solution $G$; moreover, $I_{n}-A_{0}$ is an invertible M-matrix and, hence, $(I-A_{0})^{-1}\geq 0$.

Furthermore, in view of the Perron Frobenius Theorem, there exists a unique vector $v$ such that $\mbox{\boldmath$v$}^{T}A=\mbox{\boldmath$v$}^{T}$ and $\mbox{\boldmath$v$}^{T}\mbox{\boldmath$e$}=1$, $\mbox{\boldmath$v$}>0$. If the series $\sum_{i=-1}^{\infty}iA_{i}$ is convergent we may define the vector $\mbox{\boldmath$w$}=\sum_{i=-1}^{\infty}iA_{i}\mbox{\boldmath$e$}\in\mathbb{R}^{n}$. In the study of M/G/1-type Markov chains, the drift is the scalar number $\eta=\mbox{\boldmath$v$}^{T}\mbox{\boldmath$w$}$ [14]. The sign of the drift determines the positive recurrence of the M/G/1-type Markov chain [3].

When explicitly stated, we will assume the following condition:

1. A1.

   The series $\sum_{i=-1}^{\infty}iA_{i}$ is convergent and $\eta<0$.

Under assumption [A1], the componentwise minimal nonnegative solution $G$ of equation (3) is stochastic, i.e., $G\mbox{\boldmath$e$}=\mbox{\boldmath$e$}$ ([3]). Moreover $G$ is the only stochastic solution.

In this section we reinterpret classical fixed point iterations for solving the matrix equation (3) in terms of iterative methods for solving a structured linear system.

Formally, the power series matrix equation (3) can be rewritten as the following block Toeplitz block Hessenberg linear system

The above linear system can be expressed in compact form as

where $\mbox{\boldmath$\hat{X}$}^{T}=\left[X^{T},{X^{T}}^{2},\ldots\right]^{T}$, $\tilde{P}$ is the matrix obtained from $P$ in (1) by removing its first block row and block column and $\mbox{\boldmath$E$}=\left[I_{n},0_{n},\ldots\right]^{T}$. Classical fixed point iterations for solving (3) can be interpreted as iterative methods for solving (6), based on suitable partitionings of the matrix $H$.

For instance, from the partitioning $H=M-N$, where $M=I$ and $N=\tilde{P}$, we find that the block vector $\hat{X}$ is a solution of the fixed point problem

From this equation we may generate the sequence of block vectors

such that

We may easily verify that the sequence $\{X_{k}\}_{k}$ coincides with the sequence generated by the so called natural fixed point iteration $X_{k+1}=\sum_{i=-1}^{\infty}A_{i}X_{k}^{i+1}$, $k=0,1,,\ldots$, applied to (3).

Similarly, the Jacobi partitioning, where $M=I\otimes(I_{n}-A_{0})$ and $N=M-H$, which leads to the sequence

corresponds to the traditional fixed point iteration

The anti-Gauss-Seidel partitioning, where $M$ is the block upper triangular part of $H$ and $N=M-H$, determines the $U$-based fixed point iteration

The convergence properties of these three fixed point iterations are analyzed in [3]. Among the three iterations (9) is the fastest and also the most expensive since it requires the solution of a new linear system (with multiple right hand sides) at each iteration. Moreover, it turns out that fixed-point iterations exhibit arbitrarily slow convergence for problems which are close to singularity. In particular, for positive recurrent Markov chains having a drift close to zero the convergence slows down and the number of iterations becomes arbitrarily large. In the next sections we present some new fixed point iterations which offer several advantages in terms of computational efficiency and convergence properties when compared with (9).

Recently in [7] a comparative analysis has been performed for the asymptotic convergence rates of some regular splittings of a non-singular block upper Hessenberg M-matrix of finite size. The conclusion is that the staircase splitting is faster than the anti-Gauss-Seidel splitting, that in turn is faster than the Jacobi splitting. The second result is classical, while the first one is somehow surprising since the matrix $M$ in the staircase splitting is much more sparse than the corresponding matrix in the anti-Gauss-Seidel partitioning and the splittings are not comparable. Inspired from these convergence properties, we introduce a new fixed point iteration for solving (3), based on the staircase partitioning of $H$, namely,

The splitting has attracted interest for applications in parallel computing environments [12, 10]. In principle the alternating structure of the matrix $M$ in (10) suggests several different iterative schemes.

From one hand, the odd block entries of the system $M\mbox{\boldmath$Z$}_{k+1}=N\mbox{\boldmath$\hat{X}_{k}$}+\mbox{\boldmath$E$}A_{-1}$ yield the traditional fixed point iteration. On the other hand, the even block entries lead to the implicit scheme $-A_{-1}+(I_{n}-A_{0})X_{k+1}-A_{1}X_{k+1}^{2}=\sum_{i=2}^{\infty}A_{i}X_{k}^{i+1}$ recently introduced in [4]. Differently, by looking at the structure of the matrix $M$ on the whole, we introduce the following composite two-stage iteration:

or, equivalently,

starting from an initial approximation $X_{0}$. At each step $k$, this scheme consists of a traditional fixed point iteration that computes $Y_{k}$ from $X_{k}$, followed by a cheap correction step for computing the new approximation $X_{k+1}$.

As for classical fixed point iterations, the convergence is guaranteed when $X_{0}=0$.

A similar result also holds for the case where $X_{0}$ is a stochastic matrix, assuming that [A1] holds, so that $G$ is also stochastic.

.

Propositions 1 and 2 are global convergence results. An estimate of the rate of convergence of (12) will be provided in Section 6, together with a comparison with other existing methods.

At each iteration, the scheme (12) determines two approximations which can be combined by using a relaxation technique, that is, the approximation computed at the $k$-th step takes the form of a weighted average between $Y_{k}$ and $X_{k+1}$,

In matrix terms, the resulting relaxed variant of (12) can be written as

If $\omega_{k}=0$, for $k\geq 1$, the relaxed scheme reduces to the traditional fixed point iteration (8). If $\omega_{k}=1$, for $k\geq 1$, the relaxed scheme coincides with (12). Values of $\omega_{k}$ greater than 1 can speed up the convergence of the iterative scheme.

Concerning convergence, the proof of Proposition 1 can immediately be generalized to show that the sequence $\{X_{k}\}_{k}$ defined by (13), with $X_{0}=0$, converges for any $\omega_{k}=\omega$, $k\geq 1$, such that $0\leq\omega\leq 1$. Moreover, let $\{X_{k}\}_{k\in\mathbb{N}}$ and $\{\hat{X}_{k}\}_{k\in\mathbb{N}}$ be the sequences generated by (13) for $\omega_{k}=\omega$ and $\omega_{k}=\hat{\omega}$ with $0\leq\omega\leq\hat{\omega}\leq 1$, respectively. It can be easily shown that $G\geq\hat{X}_{k}\geq X_{k}$ for any $k$ and, hence, that the iterative scheme (12) converges faster than (13) if $0\leq\omega_{k}=\omega<1$.

The convergence analysis of the modified scheme (13) for $\omega_{k}>1$ is much more involved since the choice of a relaxation parameter $\omega_{k}>1$ can destroy the monotonicity and the nonnegativity of the approximation sequence, which is at the core of the proofs of Proposition 1 and Proposition 2 . In order to maintain the convergence properties of the modified scheme we introduce the following definition.

It is worth noting that condition (14) is implicit since the construction of $X_{k+1}$ also depends on the value of $\omega_{k+1}$. By replacing $X_{k+1}$ in (14) with the expression in the right-hand side of (13) we obtain a quadratic inequality with matrix coefficients in the variable $\omega_{k+1}$. Obviously $\omega_{k}=\omega$, $0\leq\omega\leq 1$, $k\geq 1$, are eligible sequences.

The following generalization of Proposition 1 holds.

In the following section we perform a convergence analysis to estimate the convergence rate of (13) in the stationary case $\omega_{k}=\omega$, $k\geq 1$, as a function of the relaxation parameter.

Relaxation techniques are usually aimed at accelerating the convergence speed of frustratingly slow iterative solvers. Such inefficient behavior is typically exhibited when the solver is applied to a nearly singular problem. Incorporating some relaxation parameter into the iterative scheme (3) can greatly improve its convergence rate. Preliminary insights on the effectiveness of relaxation techniques applied for the solution of the fixed point problem (7) come from the classical analysis for stationary iterative solvers and are developed in Section 6.1. A precise convergence analysis is presented in Section 6.2.