0

My problem is as follows. I have Markov chain with several absorbing classes and I would like to determine which states in the state space are absorbing and which are transient. As a simple example, we could have the following sparse transition matrix, with on each row the transition probabilities of states 0,1,2,3.

import numpy as np
from scipy.sparse import csr_matrix,csgraph

P = np.array([[0.6,0.4,0.,0.],[0.3,0.7,0.,0.],[0.2,0.4,0.2,0.2],[0.,0.,0.,1.]])
P = csr_matrix(P)

Here, there are two absorbing classes: states 0-1 and state 3. State 2 is transient. Since there are classes of absorbing states, checking for probability 1 on the diagonal is not sufficient to identify the absorbing classes.

I'm wondering how I can programmatically identify the absorbing classes. Using tools from csgraph, I can find the strongly connected components using

components,labels = csgraph.connected_components(P,connection="strong")
print( components, labels )

3 [0 0 2 1]

So, I have 3 components, each with their own labels. The absorbing classes correspond certainly to one of these components, but there are also transient components. How do I find out which component is absorbing and which is transient?

Does anyone have a nice idea for this?

Forzaa
  • 1,465
  • 4
  • 15
  • 27
  • Do you mean checking whether the diagonal is `1`? – Imanol Luengo Jun 01 '16 at 12:18
  • @ImanolLuengo I mentioned that already in the question.That only works if the absorbing class consists of a single state. One could probably check whether the subgraphs of strongly connected components have row sums and column sums that add up to 1, but I fear that may be quite inefficient. – Forzaa Jun 01 '16 at 12:28
  • Without a proper definition of what does it take for an state to be *absorbing* and what does it have to satisfy to be *transient* is quite hard to help. – Imanol Luengo Jun 01 '16 at 12:35
  • @ImanolLuengo https://en.wikipedia.org/wiki/Absorbing_Markov_chain . This reference talks about a single state, but in general there can also be sets of states that once entered can never be left. This means that there is no probability flow from the set of absorbing states to the other states, whereas transient states do have positive probability of reaching absorbing states. – Forzaa Jun 01 '16 at 12:43
  • That helps thanks! If I'm not wrong, an absorbing state can be detected if it has probability to escape to other states (e.g. it has weights >0 outside the diagonal). – Imanol Luengo Jun 01 '16 at 13:04
  • Your "absorbing classes" are also known as "closed communicating classes". That might help anyone searching for background information on the problem. – Warren Weckesser Jun 01 '16 at 15:53
  • Thanks @WarrenWeckesser, now everything makes sense :) – Imanol Luengo Jun 01 '16 at 16:02

2 Answers2

2

Solving Markov chain problems where you need to rely on floating point math (np.isclose) indicates your graph representation isn't correct for the problem you're trying to solve.

Below I've presented an alternative solution using networkx (but I assume changing to csgraph wouldn't be hard) that relies solely on connectivity and working out which edges link strong components.

import networkx as nx

def absorbing_classes(g):
    internal = set()

    # Tarjan is approximately linear O(|V| + |E|)
    cpts = list(nx.strongly_connected_component_subgraphs(g))
    for sg in cpts:
        for e in sg.edges():
            internal.add(e)

    # find all the edges that aren't part of the strongly connected components
    # ~ O(E)
    transient_edges = set(g.edges()) - internal

    # find the start of the directed edge leading out from a component
    # ~ O(E)
    transient_srcs = set([ e[0] for e in transient_edges ])

    # yield everything that don't have a vertex in transient_srcs
    for sg in cpts:
        if transient_srcs - set(sg.nodes()):
            yield sg

For example, on the graph you presented:

import numpy as np
P = np.array([[0.6,0.4,0.,0.],[0.3,0.7,0.,0.],[0.2,0.4,0.2,0.2],[0.,0.,0.,1.]])

g = nx.from_numpy_matrix(P, create_using=nx.DiGraph())
for subg in absorbing_classes(g):
    print subg.nodes()

You get two absorbing classes with nodes [0, 1] in one and [3] in the other.

Andrew Walker
  • 40,984
  • 8
  • 62
  • 84
1

The following should work, noting that the total probability in an absorbing set of states equals number of states in that set:

components,labels = csgraph.connected_components(P, directed=True, connection='strong',return_labels=True)
absorbingStates = np.zeros(P.shape[0],dtype=bool)
for component in range(components):
    indices = np.where(labels==component)[0]
    n = len(indices)
    if n==1:
        probSum = P[indices,indices].sum() #np.ix_ gives an error if n==1
    else:            
        probSum = P[np.ix_(indices,indices)].sum()
    if np.isclose(probSum,n):
        absorbingStates[indices] = True

I don't think it is very efficient to constantly select submatrices from P. Better ideas are welcomed.

Forzaa
  • 1,465
  • 4
  • 15
  • 27