I try to find the stationary vector of the Markov chain by the relaxation scheme, induced to Gauss-Seidel iteration method. Hereinafter I follow Stewart W. J. (2000). Numerical Methods for Computing Stationary Distributions of Finite Irreducible Markov Chains. In Computational Probability (pp. 81-111) The method is quite simple:
- To find the infinitesimal generator $Q$ of the transition matrix.
- To decompose $Q^T=D-E-F$, where $D$ is the diagonal of $Q^T$, $-E$ is the strictly lower triangular part and $- F$ is the strictly upper triangular part of $Q^T$.
- To organize Gauss-Seidel (forward) iteration $(D - E)x^{(k+1)} = Fx^{(k)}$, $x^{(k)}$ — $k$th iteration.
The method is straightforward, takes a few of lines of code, but $x^{(k)}$ converges to 0 for the chains I tested.
In the reference above is stated, that the convergence of forward iteration is governed by the spectral radius of $(D - E)^{-1} F$. For all my chains the radius equaled 1 (no convergence, as I understand).
What's the matter?? Maybe, the method doesn't work for arbitrary chain?
If the matrix $Q^T$ is irreducible, then the Perron-Frobenius theorem states that there exists a unique vector $x$ with strictly positive elements such that $Q^Tx = 0$. Now consider trying to calculate this stationary vector. Let $Q^T = M - N$ be any splitting of $Q^T$ in which $M$ is invertible. Then the stationary iterative method
$$ x^{(k+1)} = M^{-1}Nx^{(k)},\quad k = 0,1,2,\ldots, $$
can be used to solve $Q^Tx = 0$, where $M^{-1}N = I - M^{-1}Q^T$ is called the iteration matrix. This iterative method converges to a solution of $Q^Tx = 0$ for any initial guess $x^{(0)}$ if and only if the iteration matrix $H = M^{-1}N$ is semiconvergent, that is, if and only if the limit $\lim_{k\to\infty} H^k$ exists.
Let us consider using Gauss-Seidel to solve $Q^Tx = 0$. Suppose that $Q^T$ is irreducible. Then the Gauss-Seidel splitting $Q^T = (D - E) - F$ is a regular splitting, that is, both $(D - E)^{-1}$ and $F$ are nonnegative. Now there is a theorem which says that the iteration matrix $H = (D - E)^{-1}F$ corresponding to this regular splitting is semiconvergent if all its eigenvalues other than $\lambda = 1$ belong to the open unit disk in the complex plane. Herein lies the rub. Although the irreducible matrix $Q^T$ may have $\lambda = 1$ as its only eigenvalue of maximum modulus, the Gauss-Seidel iteration matrix may not share this property, and hence may fail to converge when applied to $Q^Tx = 0$.
This situation can be remedied by defining a corresponding weighted iteration. For example, if $H$ is the Gauss-Seidel iteration matrix for $Q^T$, define the weighted iteration matrix
$$ H_\omega = (1-\omega)I + \omega H $$
Then $H_\omega$ is semiconvergent for all $\omega\in(0,1)$. Essentially, the weighting factor pulls an complex eigenvalues with modulus $1$ off the boundary of the unit disk. In general we can make the following claim:
Let $A$ be an irreducible singular M-matrix and let $A = M-N$ be any (weak) regular splitting such that $H = M^{-1}N$ is semiconvergent. Then the iterative method
$$ x^{(k+1)} = Hx^{(k)},\quad k = 0,1,2,\ldots, $$
converges to the unique (up to scaling) strictly positive solution of $Ax = 0$ for any strictly positive initial guess $x^{(0)}$.