The problem: Given an arbitrary matrix which represents the probability of transition from one state to another state in one step for a Markov chain, determine all closed sets of states and the transient state.
This small note describes an algorithm that I developed which solves this problem. The algorithm is recursive in nature.
A Matlab implementation is given below with an example run on 3 different matrices and how to call the matlab function.
This algorithm describes how to find all closed sets and the transient set (if any) given as input the \(p\) matrix which contains the initial one step finite chain Markov probability transition from one state to another. The algorithm works directly on the matrix \(p\) itself and does not require prior knowledge of which states are recurrent.
An implementation of the algorithm is provided (in Matlab) with a number of examples.
As an example of an input, let us consider the following simple Markov chain state diagram and its corresponding \(p\) matrix and then we show the output generated from this algorithm. In this diagram, the numbers on the arrows between the states (the circles) is the probability of going from the source state to the end state.
The matrix \(p\) for the above diagram is the following
The output of the algorithm below will be the following
closed sets: \(\left \{ a\right \} ,\left \{ b,e\right \} \), transient set: \(\left \{ c,d\right \} \)
Once the closed sets and the transient set are known, generating the canonical form is a simple matter of formatting the \(p\) matrix in the following form
Now the algorithm is described.
The basic idea of the algorithm is the following: For each state, we generate the path of states that can be travelled from that state. Next, we collect the states which has the same states on its path. Each set of states which has the same states on their path are one closed set. The remaining states, if any, would all have different set of states on their path. These states make up the set of transient states.
The algorithm will terminate since the length of the longest path that could possibly be travelled from any state is finite (which is the number of states in the chain).
Since there are \(n\) states, and for each state we need to examine \(n-1\) other states at most, the algorithm is of order \(O\left ( n^{2}\right ) \).
Here is a description of the algorithm
Input: A square \(n\) by \(n\) probability transition matrix (i.e. each row in the matrix sums to \(1\))
output: all sets which are closed and transient set (if it exists)
For each state \(i\) in the LHS set, scan its RHS list of states \(j\), and do the following for
each state \(j\) found
To illustrate how this algorithm works, we apply it to the example Markov chain shown above.
In step 1, we generate the LHS and the RHS lists from the input \(p\) matrix.
\begin {equation} \overbrace {\begin {Bmatrix} a\\ b\\ c\\ d\\ e \end {Bmatrix} }=\overbrace {\begin {Bmatrix} a & & & \\ b & e & & \\ b & c & d & e\\ a & c & d & e\\ b & e & & \end {Bmatrix} } \tag {1} \end {equation}
The above is send to step 2. In step 2 we start by scanning the RHS list of each state listed in the LHS set. For state \(a\) we see that there is only state \(a\) in its RHS. and from the rule in \(2.a\) we see that we do nothing and go to the next state in its RHS set. But this is the only state in this set. Hence we are done with state \(a\) and we go to state \(b\) in the LHS set.
We start again by looking at the set of states of the RHS of state \(b\), which are \(\left \{ b,e\right \} \). we start with state \(b\) which is the same as the state in the LHS, so we skip this per rule \(2.a\), and go to the next state which is \(e\). According to \(2.b\), we replace this by \(e\cup RHS\) states of \(e\), where \(RHS\) of state \(e\) is seen to be \(\left \{ b,e\right \} \) which is the last RHS set in the last entry of the LHS set as shown in (1). Hence now the RHS of state \(b\) becomes \(\left \{ b,e,\left \{ b,e\right \} \right \} \). Since there are no more states in the RHS of \(b\) we have completed this state. Now we remove duplicated entries in the above list to obtain \(\left \{ b,e\right \} \). We now go to state \(c\) in the LHS set. and do the same. We do this until we processed all states in the LHS. At the end we have the following LHS and RHS generated
\begin {equation} \overbrace {\begin {Bmatrix} a\\ b\\ c\\ d\\ e \end {Bmatrix} }=\overbrace {\begin {Bmatrix} a & & & & \\ b & e & & & \\ b & e & c & d & a\\ a & c & b & d & e\\ b & e & & & \end {Bmatrix} } \tag {2} \end {equation}
We now compare RHS sets in (1) with those in (2). Since they are not the same, we repeat step 2 again, but use the above RHS now as an input to step 2.
Again, we start by scanning the RHS list of each state in the LHS. We repeat the same thing as before, lets say we reached state \(c\) in the LHS set now. We see that its RHS is \(\left \{ b,e,c,d,a\right \} \), we start with state \(b\) and replace that with \(b\cup RHS\) states of \(b\) which is \(\left \{ b,\left \{ b,e\right \} \right \} \). We go to state \(e\), and replace that with \(e\cup RHS\) states of \(e\) which is \(\left \{ e,\left \{ b,e\right \} \right \} \). Next we do state \(c\) and replace that with \(c\cup RHS\) states of \(c\) which is \(\left \{ c,\left \{ b,e,c,d,a\right \} \right \} \), next we do state \(d\), and replace that with \(d\cup RHS\) states of \(d\) which is \(\left \{ d,\left \{ a,c,b,d,e\right \} \right \} \), and finally we do state \(a\) and replace that with \(a\cup RHS\) states of \(a\) which is \(\left \{ a,\left \{ a\right \} \right \} \), now that we completed the RHS processing of state \(c\) we remove all duplicated states in its RHS and we obtain \(\left \{ b,e,c,d,a\right \} \). Next we go to state \(d\) in the LHS set and process that similarly. We then do state \(e\) in the LHS set. Now we have completed step 2 again and the final result is the following
\begin {equation} \overbrace {\begin {Bmatrix} a\\ b\\ c\\ d\\ e \end {Bmatrix} }=\overbrace {\begin {Bmatrix} a & & & & \\ b & e & & & \\ b & e & c & d & a\\ a & c & b & d & e\\ b & e & & & \end {Bmatrix} } \tag {3} \end {equation}
We see now that the RHS in (3) is the same as the RHS in \(\left ( 2\right ) \) hence we stop and go to step 4 in the algorithm.
Now we collect all states from the LHS which has the same list in its RHS. We see that state \(a\) is its own class. And since the LHS of this class which is state \(a\) is the same as the state in the RHS which is \(a\), then this class is closed set. Next we see that states \(\left \{ b,e\right \} \) has the same states on their RHS which is \(\left \{ b,e\right \} \) and since the LHS states is the same as the RHS states, then this is a close set class that contains states \(\left \{ b,e\right \} \), next we see that state \(c\) and state \(d\) also has the same RHS set which is \(\left \{ b,e,c,d,a\right \} \) but now since \(\left \{ c,d\right \} \neq \left \{ b,e,c,d,a\right \} \), then the states \(\left \{ c,d\right \} \) make a class which is the transient states. Hence \(\left \{ c,d\right \} \) are transient and all other states are recurrent.
The following is a test run of the implementation showing in each case the matrix \(p\) and the output of the algorithm.
clear all close all nmaTestMarkov() *************************** 1.0000 0 0 0 0 0 0.2000 0.8000 0 0 0 0.7000 0.3000 0 0 0.1000 0 0.1000 0.4000 0.4000 0 0.1000 0.3000 0.2000 0.4000 found the following closed sets {1} {2,3} found the following transient set {4,5} *************************** 0.1000 0.3000 0 0 0.6000 0 0.2000 0.7000 0 0.1000 0 0.7000 0.3000 0 0 0.1000 0 0.1000 0.4000 0.4000 0 0.1000 0.3000 0.2000 0.4000 found the following closed sets {1,2,3,4,5} No transient set found *************************** 0.5000 0.5000 0 0 0 0 0 0 1.0000 0 0 0 0.3333 0 0 0.3333 0.3333 0 0 0 0 0.5000 0.5000 0 0 0 0 0 0 1.0000 0 0 0 0 1.0000 0 found the following closed sets {5,6} found the following transient set {1,2,3,4} *************************** 1.0000 0 0 0 0 0 0 1.0000 0 0 0 0 0 0 1.0000 0 0 0 0 0 0 0 0.5000 0.5000 0.2500 0.2500 0 0 0 0.5000 0.2500 0 0.2500 0 0.5000 0 found the following closed sets {1} {2} {3} found the following transient set {4,5,6} *************************** 1.0000 0 0 0 0 0 0 0 0 0 1.0000 0 0 0 0 0 0 0 0.2222 0.1111 0 0.0833 0.1111 0.1389 0.1389 0.1111 0.0833 0.1389 0.1667 0 0.7500 0 0 0 0 0 0.1111 0.1667 0 0 0.7222 0 0 0 0 0.1389 0.1667 0 0 0 0.6944 0 0 0 0.1389 0.1667 0 0 0 0 0.7222 0 0 0.1111 0.1667 0 0 0 0 0 0.7222 0 0.0833 0.1667 0 0 0 0 0 0 0.7500 found the following closed sets {1} {2} found the following transient set {3,4,5,6,7,8,9}
function nmaTestMarkov() p1=[1 0 0 0 0; 0 .2 .8 0 0; 0 .7 .3 0 0; .1 0 .1 .4 .4; 0 .1 .3 .2 .4]; testIt(p1); fprintf('\n\n'); p1=[ .1 .3 0 0 .6; 0 .2 .7 0 .1; 0 .7 .3 0 0; .1 0 .1 .4 .4; 0 .1 .3 .2 .4]; testIt(p1); fprintf('\n\n'); p1=[ .5 .5 0 0 0 0; 0 0 1 0 0 0; 1/3 0 0 1/3 1/3 0; 0 0 0 1/2 1/2 0; 0 0 0 0 0 1; 0 0 0 0 1 0]; testIt(p1); p1=[ 1 0 0 0 0 0; 0 1 0 0 0 0; 0 0 1 0 0 0; 0 0 0 0 1/2 1/2; 1/4 1/4 0 0 0 .5; 1/4 0 1/4 0 1/2 0]; %closed {1},{2},{3} testIt(p1); p1=[ 1 0 0 0 0 0 0 0 0; 0 1 0 0 0 0 0 0 0; 2/9 1/9 0 1/12 1/9 5/36 5/36 1/9 1/12; 5/36 1/6 0 3/4 0 0 0 0 0; 1/9 1/6 0 0 13/18 0 0 0 0; 5/36 1/6 0 0 0 25/36 0 0 0; 5/36 1/6 0 0 0 0 13/18 0 0; 1/9 1/6 0 0 0 0 0 13/18 0; 1/12 1/6 0 0 0 0 0 0 3/4]; %closed {1},{2},{3} testIt(p1); end %%%%%%%%%%%%%%%%%%%%%%%%%%%% % % %%%%%%%%%%%%%%%%%%%%%%%%%%%% function testIt(p) [nRow,nCol]=size(p); [closedSet,tranSet]=nmaChainMarkov(p); fprintf('***************************\n\n'); disp(p); if length(closedSet)>0 fprintf('found the following closed sets\n'); emit(closedSet); else fprintf('No closed sets found\n'); end if length(tranSet)>0 fprintf('found the following transient set\n'); emitV2(tranSet); else fprintf('No transient set found\n'); end end %%%%%%%%%%%%%% % % %%%%%%%%%%%%%% function emit(s) for i=1:length(s) c=cell2mat(s{i}(1)); fprintf('{'); for j=1:length(c) fprintf('%d',c(j)); if(j<length(c)) fprintf(','); else fprintf('}'); end end fprintf('\n'); end end %%%%%%%%%%%%%% % % %%%%%%%%%%%%%% function emitV2(s) fprintf('{'); for i=1:length(s) c=cell2mat(s{i}(1)); for j=1:length(c) fprintf('%d',c(j)); if i<length(s) fprintf(','); else if j<length(c) fprintf(','); end end end end fprintf('}\n'); end
function [cs,tr]=nmaChainMarkov(p) %function [cs,tr]=nmaChainMarkov(p) % %By Nasser Abbasi %March 8, 2008 % [nRow,nCol]=size(p); if nRow~=nCol error 'matrix must be square'; end lhs=1:nCol; rhs=cell(nRow,1)'; for i=1:nCol rhs{i}=find(p(i,:)~=0); end keepOn=1; while keepOn newRhs=updateRhs(lhs,rhs); if isequal(newRhs,rhs) keepOn=0; else rhs=newRhs; end end theClasses=findSets(lhs,newRhs); [cs,tr]=classifyClasses(theClasses); end %%%%%%%%%%%%%%%%%%%%%%%%%%%% % called after all classes are found % this function classify the class as closed or % transient %%%%%%%%%%%%%%%%%%%%%%%%%%%% function [cs,tr]=classifyClasses(c) n=length(c); cs=cell(0); tr=cell(0); csIndex=0; trIndex=0; for i=1:n if isequal( c{i}(1), c{i}(2) ) csIndex=csIndex+1; cs{csIndex}(1)=c{i}(1); cs{csIndex}(2)=c{i}(2); else trIndex=trIndex+1; tr{trIndex}(1)=c{i}(1); tr{trIndex}(2)=c{i}(2); end end end %%%%%%%%%%%%%%%%%%%%%%%%%%%% % updated the RHS lists, see paper. % %%%%%%%%%%%%%%%%%%%%%%%%%%%% function newRhs=updateRhs(lhs,rhs) for i=1:length(lhs) n=0; for j=1:length(rhs{i}) if rhs{i}(j)==lhs(i) n=n+1; tmp{i}(n)=lhs(i); else n=n+1; tmp{i}(n)=rhs{i}(j); z=rhs{i}(j); for k=1:length(rhs{z}) n=n+1; tmp{i}(n)=rhs{z}(k); end end end tmp{i}=unique(tmp{i}); end newRhs=tmp; end %%%%%%%%%%%%%%%%%%%%% % find all the states with the same rhs % %%%%%%%%%%%%%%%%%%%%%% function theClasses=findSets(lhs,rhs) n=length(lhs); skip=cell(0); nskip=0; skip{1}(1)=0; entryNumber=0; for i=1:n if isempty(find(skip{1}==i)) entryNumber=entryNumber+1; nskip=nskip+1; skip{1}(nskip)=i; tmp=cell(2,1); k=1; tmp{1}(k)=i; tmp{2}=rhs{i}; for j=1:n if j ~= i if isempty(find(skip{1}==j)) if isequal(tmp{2},rhs{j}) k=k+1; tmp{1}(k)=j; nskip=nskip+1; skip{1}(nskip)=j; end end end end if ~isempty(tmp) theClasses{entryNumber}=tmp; end end end end