In this research paper, e cient algorithms for computation of equilibrium as well as transient probability distribution of arbitrary nite state space Continuous / Discrete Time Markov Chains are proposed. The e ective idea is to solve a structured systems of linear equations e ciently. The ideas of this research paper could be utilized in solving structured system of linear equations that arise in other applications.
Markov chains provide interesting stochastic models of natural/arti cial phenomena arising in science and engineering. The existence of equilibrium behavior enables computation of equilibrium performance measures. Thus, researchers invested considerable e ort in e ciently computing the equilibrium probability distribution of Markov Chains and thus, the equilibrium performance measures.
Traditionally, computation of transient probability distribution of
Continuous Time Markov Chains (CTMCs) was considered to be a di cult open
problem. It requires computation of Matrix Exponential associated with the generator
matrix of CTMC. Even in the case of nite state space CTMCs, e cient
computation of transient probability distribution was considered to be a di cult open
problem. In [
It is well known that computation of equilibrium distribution of CTMCs
reduces to solution of linear system of equations. The approach proposed in [
This research paper is organized as follows. In Section 2, related literature is reviewed. In Section 3, algorithm for e cient computation of equilibrium probability mass function (PMF) of CTMCs is discussed. In Section 4, algorithm for e cient computation of transient PMF of CTMCs is discussed. In Section 5, we provide numerical results. The research paper concludes in Section 6. 2
It is well known that computation of equilibrium probability vector of a Homogeneous CTMC with generator Q reduces to the solution of the following linear system of equations (computation of vector in the left null space of the generator matrix Q):
Q 0: It should be noted that in the case of homogeneous Discrete Time Markov Chain (DTMC) with state transition matrix P , the computation of equilibrium probability vector reduces to the solution of following linear system of equations i.e.
P or, equivalently, (P
I)
0:
It readily follows that (P I) has the properties of a generator matrix. Thus,
in the following discussion, we consider only the computation of equilibrium
probability vector of a CTMC. In computational linear algebra, there are e
cient algorithms to solve linear system of equations [
But, from the point of view of computing , the structure of generator matrix Q
was not completely taken into account. Thus, the goal of this research paper is
to design e cient algorithm for computing , taking into account the structure
of generator matrix Q: Some related e ort dealing with Laplacian system of
equations is discussed in [
CTMC with 4 States
First consider the generator matrix of a CTMC with 4 states (i.e. 4 x 4 matrix)
partitioned using blocks of size 2 x 2 (i.e. there are 4 such 2 x 2 matrices in the
generator) [
Q =
A11 A12
A21 A22 Q~ =
A~11 A~12 :
0 A22 1A~11 = 0; 1A~12 + 2A22 = 0: 2 =
1A~12A221: where Ai;j ; i; j 2 f1; 2g, are 2 x 2 square matrices in the 4 x 4 block matrix Q.
Claim: For an irreducible positive recurrent (recurrent non-null) CTMC with generator matrix Q, it is well known that the matrices A11 and A22 are nonsingular. Hence, we have that
A21 = A22X with X = A221A21: Thus, by means of elementary column operations (i.e. multiplication by a scalar and addition/subtraction), the generator matrix Q can be converted into the block upper triangular form: Note that A~22 = A22.
Let X = [f1 f2]. Hence, column operations to arrive at Q~ (from Q) are determined by the column vectors f1; f2. Also, it follows from linear algebra, that the equilibrium probability vector is una ected by such column operations. Hence, with = [ 1 2]; we have the following system of linear equations: i.e. the boundary system of two equations and (1) (2)
The boundary system of equations corresponds to a singular matrix A~11 for an irreducible, recurrent non-null Markov chain (otherwise will be zero). Thus the boundary system of linear equations leads to a single linear equation in two variables 1 = [ 1;1 1;2]: Hence we have a linear equation of the form 1;1 + 1;2 = 0: Thus, 1;2 =
1;1 : 2 can be expressed in-terms of 1; we utilize the normalizing
e = 1; where e is a column vector of ones, to determine 1;1 and hence all the other equilibrium probabilities.
Remark 1: We realize that the inverse of a 2 x 2 matrix can be computed by inspection and essentially only requires computation of the determinant. For instance, let us consider a 2 x 2 matrix B:
Now, we generalize the above algorithm for the case, where the number of states, N = 2m, where m > 2. In such case, the generator matrix is of the following form: are all non-singular matrices, except A11 or Amm.
Proof: Refer [
Hence, as in the m = 2 case, by means of elementary column operations on the generator matrix Q, we arrive at the following block upper triangular matrix Q~:
6 0 Q~ = 666 ... The above system of linear equations is recursively solved to compute the equilibrium probability vector i.e.
for 2 j
1); j = (m 1A1j A~jj1 ~ 2A2j A~jj1 ~ j 1Aj 1j A~jj1 ~ m = 1A~1mA~m1m 2A~2mA~m1m
~ m 1Am 1mAm1m: As in the m = 2 case, using non-singular matrices on the diagonal of generator matrix, Q, it can be converted to a block upper-triangular matrix. Since inverse of 2 x 2 matrices on the diagonal of Q can be computed e ciently, an e cient algorithm exists for such a computational problem.
De nition: nite memory recursion of order L for the equilibrium
probability vector is of the following form
(L + 1) = (1)W1 + (2)W2 +
+ (L)WL;
where Wi; i = 1; : : : ; L, are recursion matrices. We call such a recursion forward
nite memory recursion [
Remark 2: Since A~jj are all non-singular 2 x 2 matrices, their inverse can be e ciently computed.
Computational Complexity We now compute the computational complexity of the algorithm to determine equilibrium probabilities in most general case. { Conversion of Q into a Block Upper Triangular Matrix: This task is accomplished using an approach similar to Gaussian elimination using the 2 x 2 non-singular matrices. It requires (i) (M 1)2 + (M 2)2 + + 22 + 1 =
M (M 1)(2M 6 1) ; 2 x 2 matrix multiplication,
M (M 1)(2M 6 1) ; 2 x 2 matrix additions and (M 1) inversions of 2 x 2 matrices. (ii) 2 x 2 Matrix inversions needed in Finite Memory Recursion: We require (m 2) inversions of 2 x 2 matrices (since Am~ m is already inverted in step (i)). Each 2 x 2 matrix inversion requires (a) Determinant computation i.e. 2 multiplications and one subtraction; (b) Division of elements by determinant i.e. 4 divisions. (iii) 2 x 2 Matrix Multiplications: In the above system of linear equations, we require 1 + 2 + 3 + + (m 1) = 2 multiplications of two 2 x 2 matrices. Each such multiplication (of two 2 x 2 matrices) requires at most 8 multiplications and 4 additions. To reduce the complexity, we can utilize Strassen's multiplication algorithm. The number of multiplications required to compute the product of 2 x 2 matrices reduces to 7 (using Strassen's algorithm). But with Strassen's algorithm we require 18 additions/subtractions. (iv) Vector Matrix Multiplications: We require m(m
1) 1 + 2 + 3 + + (m 1) = m(m
1) 2 multiplications of 1 x 2 and 2 x 2 matrices ( i.e. multiplication of row vector and 2 x 2 matrix). { Normalization requires (2m 1) additions and one division to obtain 1;1: { Finally, we require (2m 1) multiplications to determine the equilibrium probabilities using 11: Thus, total number of 2 x 2 matrix inversions, 2 x 2 matrix multiplications and 2 x 2 matrix additions and other arithmetical operations required can be readily determined.
Note that in the computer implementation of any algorithm, divisions are \more expensive" than multiplications. Hence, in the above inversion of 2 x 2 matrices, we rst compute the inverse of determinant and convert the associated divisions to multiplications.
Note also that in the case of QBD processes, Latouche et.al utilize Gaussian
Elimination method to compute the equilibrium probability distribution [
Remark 3: Suppose, the number of states of the CTMC is an Odd number. Then, we consider the 3 x 3 boundary system of linear equations (i.e. the initial probability vector is of dimension 3) and utilize the above idea to compute the equilibrium probability vector e ciently.
Remark 4: It should be noted that using similar idea (as discussed above), by
means of elementary column operations, the generator matrix can be converted
into a block lower triangular matrix. In this case, the boundary system of linear
equations is at the trailing boundary. By solving for the last probability vector,
we recursively compute the equilibrium probabilities. We call such a recursion
as the backward nite memory recursion [
It is well known that the transient PMF of a homogeneous CTMC (generator matrix doesn't depend on time t, unlike non-homogeneous CTMC) satis es the following vector matrix di erential equation i.e.
d dt
(t) = (t)Q: s~(s)
~(0) = ~(s)Q: ~(s) = ~(0)[Q sI] 1: Taking Laplace-Transform on both sides, we have that Equivalently, we have the following expression for ~(s): These constitute a structured system of linear equations. Our goal is to e ciently solve such system of equations for (s) and compute the inverse Laplace transform of (s) to arrive at the time dependent (transient) PMF of CTMC.
Lemma 2: For a positive recurrent (recurrent non-null) CTMC, in the Region of Convergence (ROC) (of Laplace transform) the sub-matrices on the diagonal of (Q sI) are non-singular.
Proof : As in the equilibrium case, such matrices are all strictly diagonally dominant and hence are all non-singular. tu
Thus, as in the case of computation of equilibrium PMF, (s) can be determined e ciently for the values of s lying in the Region of Convergence. Detailed duplication of equations is avoided for brevity.
In this research paper, e cient algorithms for computing the equilibrium and transient probability distribution of an arbitrary nite state space CTMC are discussed. The algorithms e ectively solve a structured system of linear equations e ciently. It might be interesting to apply the proposed method to solving system of linear equations, where the coe cient matrix has the structure similar to a generator matrix (i.e. structured diagonally dominant matrix), for instance, the generator matrix has Toeplitz structure in some applications. This idea, however, is left for future research. 12. Stothers, Andrew James. "On the complexity of matrix multiplication." (2010). 13. Strassen, Volker. "Gaussian elimination is not optimal." Numerische mathematik 134), 354-356 (1969). 14. Strassen, Volker. "The asymptotic spectrum of tensors and the exponent of matrix multiplication." In 27th Annual Symposium on Foundations of Computer Science (sfcs 1986), pp. 49-54. IEEE, (1986). 15. Williams, Virginia Vassilevska. "Multiplying matrices faster than CoppersmithWinograd extended abstract." In Proceedings of the 2012 ACM Symposium on Theory of Computing (STOC'12). ACM, New York. (2012). 16. Zhang, Ji, and Edward J. Coyle. "Transient analysis of quasi-birth-death processes." Stochastic Models 5(3), 459-496 (1989).