A Numerical Algorithm on the Computation of the Stationary Distribution of a Discrete Time Homogenous Finite Markov Chain

Abstract

The transition matrix, which characterizes a discrete time homogeneous Markov chain, is a stochastic matrix. A stochastic matrix is a special nonnegative matrix with each row summing up to 1. In this paper, we focus on the computation of the stationary distribution of a transition matrix from the viewpoint of the Perron vector of a nonnegative matrix, based on which an algorithm for the stationary distribution is proposed. The algorithm can also be used to compute the Perron root and the corresponding Perron vector of any nonnegative irreducible matrix. Furthermore, a numerical example is given to demonstrate the validity of the algorithm.

1. Introduction and Preliminaries

Throughout this paper, the following notations and definitions are used. A matrix 𝐴=(𝑎_𝑖𝑗)∈ℝ^𝑚×𝑛 is called nonnegative (positive), if all 𝑎_𝑖𝑗≥0 (𝑎_𝑖𝑗>0), denoted by 𝐴≥0 (𝐴>0). Similarly, a vector 𝑥=(𝑥₁,…,𝑥_𝑛)^𝑇 is called nonnegative (positive) and denoted by 𝑥≥0 (𝑥>0), if all 𝑥_𝑖≥0 (𝑥_𝑖>0). Let 𝐵=(𝑏_𝑖𝑗)∈ℝ^𝑚×𝑛, we denote 𝐴≥𝐵 (𝐴>𝐵), if 𝐴−𝐵≥0 (>0), that is, 𝑎_𝑖𝑗≥𝑏_𝑖𝑗 (𝑎_𝑖𝑗>𝑏_𝑖𝑗) for all 𝑖=1,2,…,𝑚,𝑗=1,2,…,𝑛.

For a square matrix 𝐴∈ℝ^𝑛×𝑛 with eigenvalues 𝜆₁,…,𝜆_𝑛,𝜌(𝐴)=max{|𝜆_𝑗|} is called the spectral radius of 𝐴. If 𝐴≥0 is irreducible, there exists a unique eigenvector 𝑥=(𝑥₁,…,𝑥_𝑛)^𝑇>0 such that 𝐴𝑥=𝜌(𝐴)𝑥 and ‖𝑥‖₁=|𝑥₁|+⋯+|𝑥_𝑛|=1. In this case, we say that 𝜌(𝐴) is the Perron root of 𝐴 and 𝑥 is the Perron vector [1].

We consider a discrete-time Markov chain 𝑋={𝑋^(𝑛)∶𝑛=0,1,…} with a finite state space 𝒮={𝑖₁,…,𝑖_𝑛}. Among ergodic processes, homogeneous Markov chains with finite state space are particularly interesting examples. Such processes satisfy the Markov property, which states that their future behavior, conditional to the past and present, depends only on the present. Precisely, for all 𝑡∈ℝ₊, ℎ>0, and for all sequences 0≤𝑡₁≤⋯≤𝑡_𝑟=𝑡,𝑖₁,…,𝑖_𝑟∈𝒮 and 𝑖_𝑗∈𝒮, 𝑃𝑋^(𝑡+ℎ)=𝑖_𝑗∣𝑋^(𝑡)=𝑖_𝑟,𝑋^{(𝑡𝑟−1)}=𝑖_𝑟−1,…,𝑋^(𝑡1)=𝑖₁=𝑃𝑋^(ℎ)=𝑖_𝑗∣𝑋⁽⁰⁾=𝑖_𝑟. (1.1)

The behavior of such a process is characterized by an 𝑛×𝑛 matrix 𝑀 called the transition matrix [2].

Its stationary distribution 𝜋, which is also its asymptotic distribution, is a vector satisfying the following. 𝜋^𝑇𝑀=𝜋^𝑇,^𝑛_𝑗=1𝜋_𝑗=1, (1.2) that is, 𝜋^𝑇𝑀=𝜋^𝑇,𝜋>0,𝜋^𝑇𝑒=1, (1.3) where 𝑒 is the column vector of all ones.

It has been established that it is possible to represent all possible uses of a software system as a Markov chain [3–5]. This model is called a Markov chain usage model. In a usage model, states of use (such as state “Document Loaded” in a model representing use of a word processing system) are represented by states in the Markov chain. Transitions between states of use (such as moving from state “Document Loaded” to “No Document Loaded” when the user closes a document in a word processing system) are represented by state transitions between the appropriate states in the Markov chain. Transitions between states of use have associated probabilities which represent the probability of making each transition. A usage model may be created based on information taken from functional specifications, usage specifications, and test objectives.

Considering the problem of software reliability, we represent a software system 𝑆_𝑓 with 𝑛 states of use {𝑠₁,…,𝑠_𝑛} by a homogeneous discrete Markov chain {𝑋^(𝑛)∶𝑛=0,1,…} (the corresponding transition matrix is 𝑀). We denote the initial state probability distribution 𝜋⁽⁰⁾=(𝜋₁⁽⁰⁾,…,𝜋_𝑛⁽⁰⁾)^𝑇, where 𝜋_𝑖⁽⁰⁾=𝑃(𝑋⁽⁰⁾=𝑠_𝑖). Then (𝜋^(𝑘))^𝑇=(𝜋^(𝑘−1))^𝑇𝑀, where 𝜋^(𝑘) stands for the state probability distribution at time 𝑘. Let 𝜇_𝑖 (𝑖=1,…,𝑛) be the probability when the software fails at state 𝑠_𝑖. The reliability of 𝑆_𝑓 at time 𝑘 can be defined as 𝑅^(𝑘)=1−∑^𝑛_𝑖=1𝜇_𝑖𝜋_𝑖^(𝑘). After a long time running, the state distribution of system 𝑆_𝑓 will tend to the stationary distribution 𝜋=(𝜋₁,…,𝜋_𝑛)^𝑇. Then, the terminating reliability 𝑅=1−∑^𝑛_𝑖=1𝜇_𝑖𝜋_𝑖, with which we can evaluate the quality of a software system. By decreasing the 𝜇_𝑖 of state 𝑠_𝑖 with the largest probability 𝜋_𝑖 in the stationary distribution 𝜋, we can also enhance the reliability of 𝑆_𝑓 efficiently with limited resources.

A nonnegative matrix 𝐴=(𝑎_𝑖𝑗)∈ℝ^𝑛×𝑛 is called a row-stochastic matrix (or a stochastic matrix for short) if ∑^𝑛_𝑗=1𝑎_𝑖𝑗=1 for all 𝑖=1,…,𝑛, that is, 𝐴𝑒=𝑒.

From the well-known Perron-Frobenius theorem, it can be easily deduced that the Perron root of a stochastic matrix 𝐴 equals 1.

Obviously, the transition matrix 𝑀 of a discrete-time homogeneous Markov chain is a stochastic matrix. From (1.3), we have 𝑀^𝑇𝜋=𝜋. That is to say, the stationary distribution 𝜋 is also an eigenvector of 𝑀^𝑇 associated to 1. Since 𝑀^𝑇 and 𝑀 have the same eigenvalues, 𝜋 is the Perron root of 𝑀^𝑇, that is, the solution to 𝑀^𝑇𝜋=𝜋. As for the computational aspects of 𝜋, many approaches have been presented (e.g., see [6–10]) based on the Gaussian elimination, direct projection and so on. In this paper, from the viewpoint of the Perron root which has not been discussed, we propose an algorithm for the stationary distribution 𝜋 considering that the computation of 𝜋 is equivalent to the computation of the Perron vector of 𝑀, which not only can compute the stationary distribution, but also could be used to compute the Perron root and the corresponding Perron vector of any nonnegative irreducible matrix (noting that the stationary distribution is the Perron vector of the transition matrix, which is a special nonnegative matrix).

This paper is organized as follows. In the next section, we propose some lemmas and preliminary results. In Section 3, we prove the convergent theorem and give some facts. In Section 4, we propose an algorithm for the stationary distribution together with a demonstrating numerical example.

2. Some Lemmas

In this section, we present some lemmas which will be used in the proof of the main results. The following facts can be found in [1, 11, 12].

Definition 2.1 (see [1]). Let 𝐴∈ℝ^𝑛×𝑛 be a nonnegative matrix. If 𝐴^𝑚>0 for some integer 𝑚≥1, one says that 𝐴 is primitive.

It is known that any primitive matrix must be irreducible [12]. We will use the following important facts which can be found in [1, 12].

Theorem A (see [1]). Let 𝐴=𝐴_𝑛×𝑛≥0, then 𝐴 is irreducible if and only if (𝐼+𝐴)^𝑛−1>0, where 𝐼 is the unit matrix.

Theorem B (Perron-Frobenius (see [1])). Let 𝐴=𝐴_𝑛×𝑛≥0 be irreducible. Then, (a)𝜌(𝐴)>0, (b)𝜌(𝐴) is an eigenvalue of 𝐴, (c)There exists a vector 𝑥>0 such that 𝐴𝑥=𝜌(𝐴)𝑥,(d)𝜌(𝐴) is a simple eigen value of 𝐴.

This theorem guarantees the eigenspace of 𝜌(𝐴) is one-dimensional. That is, 𝐴𝑦=𝜌(𝐴)𝑦 implies 𝑦=𝑘𝑥. And there exists an unique positive vector 𝑥>0 whose components sum to 1 such that 𝐴𝑥=𝜌(𝐴)𝑥. This 𝑥 is called the Perron vector [1].

For the Perron root of nonnegative matrices, many algorithms and bounds estimations have been proposed (see in [13, 14]). In this paper, we will describe the Perron root by using the following Collatz-Wielandt functions [11, 12].

Definition 2.2 (see [11]). Let 𝐴=(𝑎_𝑖𝑗)_𝑛×𝑛 be nonnegative, define 𝑓_𝐴(𝑥)()=min𝐴𝑥_𝑗𝑥_𝑗,𝑔_𝐴(𝑥)()=max𝐴𝑥_𝑗𝑥_𝑗, (2.1) for any positive vector 𝑥=(𝑥₁,…,𝑥_𝑛)^𝑇>0.

𝑓_𝐴(𝑥) and 𝑔_𝐴(𝑥) are both continuous at any 𝑥>0.

Lemma 2.3. Let 𝐴=𝐴_𝑛×𝑛 be nonnegative and irreducible. Then, for any 𝑥>0, 𝑓_𝐴(𝑥) and 𝑔_𝐴(𝑥) satisfy the following: (1)𝐴𝑥≥𝑓_𝐴(𝑥)𝑥,𝐴𝑥≤𝑔_𝐴(𝑥)𝑥,and𝑓_𝐴(𝑥)≤‖𝐴‖₁, (2)𝑓_𝐴(𝑡𝑥)=𝑓_𝐴(𝑥),𝑔_𝐴(𝑡𝑥)=𝑔_𝐴(𝑥)(𝑡>0), (3)𝐴𝑥>𝑡𝑥 gives 𝑓_𝐴(𝑥)>𝑡; 𝐴𝑥<𝑡𝑥 gives 𝑔_𝐴(𝑥)<𝑡,(4)If 𝐵≥0 is irreducible and 𝐴𝐵=𝐵𝐴, let 𝑦=𝐵𝑥, then 𝑓_𝐴(𝑦)≥𝑓_𝐴(𝑥) and 𝑔_𝐴(𝑦)≤𝑔_𝐴(𝑥).

Proof. (1)–(3) are clearly true (see [10]). For (4), by 𝐴𝑥≥𝑓_𝐴(𝑥)𝑥, it follows that 𝐴𝑦=𝐵(𝐴𝑥)≥𝐵𝑓_𝐴(𝑥)𝑥=𝑓_𝐴(𝑥)𝑦. This gives 𝑓_𝐴(𝑦)≥𝑓_𝐴(𝑥). Similarly, 𝑔_𝐴(𝑦)≤𝑔_𝐴(𝑥).

Lemma 2.4 (see [1]). If 𝐵≥0 is primitive (𝐵^𝑚>0 for some 𝑚≥1), then lim_𝑘→∞𝜌(𝐵)⁻¹𝐵^𝑘=𝐿>0,𝐵𝐿=𝜌(𝐵)𝐿. (2.2)

3. Main Results

In this section, we will present the main results.

Theorem 3.1. If 𝐴=𝐴_𝑛×𝑛≥0 is irreducible and 𝐵=𝐵_𝑛×𝑛≥0 is primitive such that 𝐴𝐵=𝐵𝐴. Let 𝑥⁽⁰⁾=(𝑎₁,…,𝑎_𝑛)^𝑇>0. Define for 𝑘=1,2,…, 𝑦^(𝑘)=𝐵𝑥^(𝑘−1),𝑥^(𝑘)=𝑦^(𝑘)𝑦^(𝑘)₁. (3.1) Then, (a)lim_𝑘→∞𝑥^(𝑘)=𝑥>0, and 𝐴𝑥=𝜌(𝐴)𝑥 with ‖𝑥‖₁=1, (b)lim_𝑘→∞𝑓_𝐴(𝑥^(𝑘))=lim_𝑘→∞𝑔_𝐴(𝑥^(𝑘))=𝜌(𝐴),(c)𝑓_𝐴(𝑥⁽⁰⁾)≤𝑓_𝐴(𝑥⁽¹⁾)≤⋯≤𝜌(𝐴)≤⋯≤𝑔_𝐴(𝑥⁽¹⁾)≤𝑔_𝐴(𝑥⁽⁰⁾).

Proof. By (3.1), we can write 𝑥^(𝑘)=𝐵^𝑘𝑥⁽⁰⁾/𝑏_𝑘 (for some 𝑏_𝑘>0). This means for 𝑘=1,2,…𝐵^𝑘𝑥⁽⁰⁾𝑏_𝑘1=𝑥^(𝑘)₁=1. (3.2)
By Lemma 2.4, lim_𝑘→∞𝐵𝜌(𝐵)^𝑘=𝐿>0,𝐵𝐿=𝜌(𝐵)𝐿. (3.3) Equations (3.2) and (3.3) imply that lim_𝑘→∞𝑏_𝑘𝜌(𝐵)^𝑘=𝐿𝑥⁽⁰⁾₁,lim_𝑘→∞𝑥^(𝑘)=lim_𝑘→∞𝐵^𝑘𝑥⁽⁰⁾𝑏_𝑘=𝐿𝑥⁽⁰⁾𝐿𝑥⁽⁰⁾₁. (3.4) By putting 𝑥=𝐿𝑥⁽⁰⁾/‖𝐿𝑥⁽⁰⁾‖₁, it is clear that 𝑥>0 (with ‖𝑥‖₁=1) and 𝐵𝑥=𝐵𝐿𝑥⁽⁰⁾𝐿𝑥⁽⁰⁾₁=𝜌(𝐵)𝐿𝑥⁽⁰⁾𝐿𝑥⁽⁰⁾₁=𝜌(𝐵)𝑥. (3.5) Since 𝐴𝐵=𝐵𝐴, we get 𝐵(𝐴𝑥)=𝜌(𝐵)𝐴𝑥. The Perron-Frobenius theorem (Theorem B) guarantees that 𝜌(𝐵) is a simple eigenvalue of 𝐵. So, 𝐵(𝐴𝑥)=𝜌(𝐵)𝐴𝑥 gives that 𝐴𝑥=𝜆𝑥(𝑥>0), which implies that 𝜆=𝜌(𝐴) and 𝐴𝑥=𝜌(𝐴)𝑥. On the other hand, by Definition 2.2, 𝐴𝑥=𝜌(𝐴)𝑥(𝑥>0) gives 𝑓_𝐴(𝑥)=𝜌(𝐴)=𝑔_𝐴(𝑥). By lim_𝑘→∞𝑥^(𝑘)=𝑥(𝑥>0), we conclude that lim_𝑘→∞𝑓_𝐴𝑥^(𝑘)=𝑓_𝐴(𝑥)=𝜌(𝐴),lim_𝑘→∞𝑔_𝐴𝑥^(𝑘)=𝑔_𝐴(𝑥)=𝜌(𝐴). (3.6) By Lemma 2.3 and (3.1), we have for 𝑘=1,2,…𝑓_𝐴𝑥^(𝑘−1)≤𝑓_𝐴𝑥^(𝑘),𝑔_𝐴𝑥^(𝑘)≤𝑔_𝐴𝑥^(𝑘−1). (3.7) So, {𝑓_𝐴(𝑥^(𝑘))} and {𝑔_𝐴(𝑥^(𝑘))} are both monotonic convergent sequences. This proves (c), completing the proof.

Remark 3.2. From the proof, we know 𝐴𝐵𝑥=𝜌(𝐴)𝜌(𝐵)𝑥 (𝑥>0), and 𝜌(𝐴𝐵)=𝜌(𝐴)𝜌(𝐵).
For an 𝑛×𝑛 irreducible matrix 𝐴≥0, since (𝑏𝐼+𝐴)^𝑛−1>0(𝑏>0), 𝐵=(𝑏𝐼+𝐴)^𝑚 are primitive for 𝑚=1,2,…. Clearly, 𝐵𝐴=𝐴𝐵, we have the following.

Corollary 3.3. If 𝐴=𝐴_𝑛×𝑛≥0 is irreducible and let 𝐵=(𝑏𝐼+𝐴)^𝑚 (for fixed 𝑚≥1 and 𝑏>0). Let 𝑥⁽⁰⁾=(𝑎₁,…,𝑎_𝑛)^𝑇>0. For all 𝑘=1,2,…, define 𝑦^(𝑘)=𝐵𝑥^(𝑘),𝑥^(𝑘)=𝑦^(𝑘)𝑦^(𝑘)₁. (3.8) Then, (a)lim_𝑘→∞𝑥^(𝑘)=𝑥>0,𝐴𝑥=𝜌(𝐴)𝑥,(‖𝑥‖₁=1), (b)lim_𝑘→∞𝑓_𝐴(𝑥^(𝑘))=lim_𝑘→∞𝑔_𝐴(𝑥^(𝑘))=𝜌(𝐴), (c)𝑓_𝐴(𝑥⁽⁰⁾)≤𝑓_𝐴(𝑥⁽¹⁾)≤⋯≤𝜌(𝐴)≤⋯≤𝑔_𝐴(𝑥⁽¹⁾)≤𝑔_𝐴(𝑥⁽⁰⁾).

For a positive matrix 𝐴>0, all the matrices 𝐵=𝐴^𝑚(𝑚=1,2,…) are primitive. The following is obvious.

Corollary 3.4. If 𝐴=𝐴_𝑛×𝑛>0 and 𝐵=𝐴^𝑚 (for 𝑚≥1). Let 𝑥⁽⁰⁾=(𝑎₁,…,𝑎_𝑛)^𝑇 be a positive vector. Define 𝑦^(𝑘)=𝐵𝑥^(𝑘),𝑥^(𝑘)=𝑦^(𝑘)𝑦^(𝑘)₁∀𝑘=1,2,…. (3.9) Then, (a)lim_𝑘→∞𝑥^(𝑘)=𝑥>0,𝐴𝑥=𝜌(𝐴)𝑥,(‖𝑥‖₁=1), (b)lim_𝑘→∞𝑓_𝐴(𝑥^(𝑘))=lim_𝑘→∞𝑔_𝐴(𝑥^(𝑘))=𝜌(𝐴), (c)𝑓_𝐴(𝑥⁽⁰⁾)≤𝑓_𝐴(𝑥⁽¹⁾)≤⋯≤𝜌(𝐴)≤⋯≤𝑔_𝐴(𝑥⁽¹⁾)≤𝑔_𝐴(𝑥⁽⁰⁾).

By (3.1), let 𝜆^(𝑘)=(1/2)(𝑓_𝐴(𝑥^(𝑘))+𝑔_𝐴(𝑥^(𝑘))) for (𝑘=1,2,…), one has the following.

Corollary 3.5. If 𝑔_𝐴(𝑥^(𝑘))−𝑓_𝐴(𝑥^(𝑘))<𝜀(𝜀>0), then |𝜌(𝐴)−𝜆^(𝑘)|<(1/2)𝜀.

Proof. From Theorem 3.1, it follows that 𝜌(𝐴)−𝜆^(𝑘)≤12𝜌(𝐴)−𝑓_𝐴𝑥^(𝑘)+12𝜌(𝐴)−𝑔_𝐴𝑥^(𝑘)=12𝜌(𝐴)−𝑓_𝐴𝑥^(𝑘)+12𝑔_𝐴𝑥^(𝑘)(𝐴)=1−𝜌2𝑔_𝐴𝑥^(𝑘)−𝑓_𝐴𝑥^(𝑘)<12𝜀. (3.10)

If 𝐴≥0 is irreducible, it is obvious that 𝐵=𝑏𝐼+𝐴(𝑏>0) is primitive, and 𝜌(𝐴)=𝜌(𝐵)−𝑏. So, we have the following.

Corollary 3.6. If 𝐴≥0 is irreducible, for any 𝑏>0, let 𝐵=𝑏𝐼+𝐴. Using the sequences {𝑓_𝐵(𝑥^(𝑘))} and {𝑔_𝐵(𝑥^(𝑘))}. Then, lim_𝑘→∞𝑓_𝐵𝑥^(𝑘)=lim_𝑘→∞𝑔_𝐵𝑥^(𝑘)=𝜌(𝐵);lim_𝑘→∞𝑥^(𝑘)‖𝑥‖=𝑥>0,₁=1;𝜌(𝐴)=𝜌(𝐵)−𝑏. (3.11)

4. An Algorithm and a Numerical Example

In this section, we propose a numerical algorithm to compute the stationary distribution of a discrete time homogeneous finite Markov chain.

Algorithm 4.1 (to compute the stationary distribution 𝜋). Step 1. Giving a transition matrix 𝑀 of a discrete time homogeneous finite Markov chain, a calculation precision 𝜀>0. Choosing parameters: a positive real number 𝑏>0 and an integer 𝑚. Setting the initial iterative vector 𝜋⁽⁰⁾=(1,1,…,1)^𝑇, 𝐵=(𝑏𝐼+𝑀^𝑇)^𝑚, 𝑘=1.
Step 2. Computing 𝜋^(𝑘) from 𝜋^(𝑘−1): 𝑦^(𝑘)=𝐵𝜋^(𝑘−1),𝜋^(𝑘)=𝑦^(𝑘)‖𝑦^(𝑘)‖₁(∀𝑘=1,2,…). (4.1)
Step 3. Compute 𝑓_𝑀𝑇(𝜋^(𝑘)) and 𝑔_𝑀𝑇(𝜋^(𝑘)): 𝑓_𝑀𝑇𝜋^(𝑘)=min_{1≤𝑖≤𝑛}𝑀^𝑇𝜋^(𝑘)_𝑖𝜋_𝑖^(𝑘),𝑔_𝑀𝑇𝜋^(𝑘)=max_{1≤𝑖≤𝑛}𝑀^𝑇𝜋^(𝑘)_𝑖𝜋_𝑖^(𝑘). (4.2)
Step 4. If 𝑔_𝑀𝑇(𝜋^(𝑘))−𝑓_𝑀𝑇(𝜋^(𝑘))<𝜀, go to Step 5. Otherwise setting 𝑘∶=𝑘+1, go back to Step 2.
Step 5. Let 𝜆=(1/2)(𝑓_𝑀𝑇(𝜋^(𝑘))+𝑔_𝑀𝑇(𝜋^(𝑘))). Then 𝜆 is the approximation of the Perron root of 𝑀^𝑇, and the corresponding 𝜋^(𝑘) is the approximation of the stationary distribution of 𝑀.

Remark 4.2. From Theorem 3.1, the convergence of Algorithm 4.1 is obvious.

We next give a numerical example.

Example 4.3. For a given finite Markov Chain, with the corresponding transition matrix as the following: 𝑀=0.200.30.100.400.10.200.50.20.5000.100.400.30.20.20.10.20.30.1000.40.20.200.30.500, (4.3) finding its approximating stationary distribution with calculation precision 𝜀=10⁻⁶.

By choosing the initial iterative vector 𝑥⁽⁰⁾=(1,1,1,1,1,1)^𝑇, parameters 𝑚=3,𝑏=1, that is, 𝐵=(𝐼+𝑀^𝑇)³, and applying Algorithm 4.1, the approximating Perron root and Perron vector are obtained after 12 iterations: 𝜋⁽¹²⁾=(0.2115999,0.0758786,0.18776110,0.1956869,0.0958466,0.2332268)^𝑇,𝜆⁽¹²⁾=1.0000001. (4.4) The iteration results are listed in Table 1.

Table 1 

Iteration results of Example 4.3 by Algorithm 4.1.

Acknowledgment

The project was supported by the National Natural Science Foundation of China (Grant no. 60831001).