Directed probability graph - algorithm to reduce cycles?

Question

Consider a directed graph which is traversed from first node 1 to some final nodes (which have no more outgoing edges). Each edge in the graph has a probability associated with it. Summing up the probabilities to take each possible path towards all possible final nodes returns 1. (Which means, we are guaranteed to arrive at one of the final nodes eventually.)

The problem would be simple if loops in the graph would not exist. Unfortunately rather convoluted loops can arise in the graph, which can be traversed an infinite amount of times (probability decreases multiplicatively with each loop traversal, obviously).

Is there a general algorithm to find the probabilities to arrive at each of the final nodes?

A particularly nasty example:

We can represent the edges as a matrix (probability to go from row (node) x to row (node) y is in the entry (x,y))

{{0, 1/2, 0, 1/14, 1/14, 0, 5/14}, 
 {0, 0, 1/9, 1/2, 0, 7/18, 0}, 
 {1/8, 7/16, 0, 3/16, 1/8, 0, 1/8}, 
 {0, 1, 0, 0, 0, 0, 0}, 
 {0, 0, 0, 0, 0, 0, 0}, 
 {0, 0, 0, 0, 0, 0, 0}, 
 {0, 0, 0, 0, 0, 0, 0}}

Or as a directed graph:

The starting node 1 is blue, the final nodes 5,6,7 are green. All edges are labelled by the probability to traverse them when starting from the node where they originate.

This has eight different paths from starting node 1 to the final nodes:

{{1/14, {1, 5}}, {5/14, {1, 7}}, {7/36, {1, 2, 6}}, 
 {1/144, {1, 2, 3, 5}}, {1/144, {1, 2, 3, 7}}, 
 {1/36, {1, 4, 2, 6}}, {1/1008, {1, 4, 2, 3, 5}}, {1/1008, {1, 4, 2, 3, 7}}}

(The notation for each path is {probability,sequence of nodes visited})

And there are five distinct loops:

{{1/144, {2, 3, 1}}, {7/144, {3, 2}}, {1/2, {4, 2}}, 
{1/48, {3, 4, 2}}, {1/1008, {4, 2, 3, 1}}}

(Notation for loops is {probability to traverse loop once,sequence of nodes visited}).

If only these cycles could be resolved to obtain an effectively tree like graph, the problem would be solved.

Any hint on how to tackle this?

I'm familiar with Java, C++ and C, so suggestions in these languages are preferred.

Answer 1

I'm not expert in the area of Markov chains, and although I think it's likely that algorithms are known for the kind of problem you present, I'm having difficulty finding them.

If no help comes from that direction, then you can consider rolling your own. I see at least two different approaches here:

1. Simulation.

Examine how the state of the system evolves over time by starting with the system in state 1 at 100% probability, and performing many iterations in which you apply your transition probabilities to compute the probabilities of the state obtained after taking a step. If at least one final ("absorbing") node can be reached (at non-zero probability) from every node, then over enough steps, the probability that the system is in anything other than a final state will decrease asymptotically toward zero. You can estimate the probability that the system ends in final state S as the probability that it is in state S after n steps, with an upper bound on the error in that estimate given by the probability that the system is in a non-final state after n steps.

As a practical matter, this is the same is computing Trⁿ, where Tr is your transition probability matrix, augmented with self-edges at 100% probability for all the final states.

2. Exact computation.

Consider a graph, G, such as you describe. Given two vertices i and f, such that there is at least one path from i to f, and f has no outgoing edges other than self-edges, we can partition the paths from i to f into classes characterized by the number of times they revisit i prior to reaching f. There may be an infinite number of such classes, which I will designate C_if(n), where n represents the number of times the paths in C_if(n) revisit node i. In particular, C_ii(0) contains all the simple loops in G that contain i (clarification: as well as other paths).

The total probability of ending at node f given that the system traverses graph G starting at node i is given by

Pr(f|i, G) = Pr(C_if(0)|G) + Pr(C_if(1)|G) + Pr(C_if(2)|G) ...

Now observe that if n > 0 then each path in C_if(n) has the form of a union of two paths c and t, where c belongs to C_ii(n-1) and t belongs to C_if(0). That is, c is a path that starts at node i and ends at node i, passing through i n-1 times between, and t is a path from i to f that does not pass through i again. We can use that to rewrite our probability formula:

Pr(f|i,G) = Pr(C_if(0)|G) + Pr(C_ii(0)|G) * Pr(C_if(0)|G) + Pr(C_ii(1)|G) * Pr(C_if(0)|G) + ...

But note that every path in C_ii(n) is a composition of n+1 paths belonging to C_ii(0). It follows that Pr(C_ii(n)|G) = Pr(C_ii(0)|G)ⁿ⁺¹, so we get

Pr(f|i) = Pr(C_if(0)|G) + Pr(C_ii(0)|G) * Pr(C_if(0)|G) + Pr(C_ii(0)|G)² * Pr(C_if(0)|G) + ...

And now, a little algebra gives us

Pr(f|i,G) - Pr(C_if(0)|G) = Pr(C_ii(0)|G) * Pr(f|i,G)

, which we can solve for Pr(f|i,G) to get

Pr(f|i,G) = Pr(C_if(0)|G) / (1 - Pr(C_ii(0)|G))

We've thus reduced the problem to one in terms of paths that do not return to the starting node, except possibly as their end node. These do not preclude paths that have loops that don't include the starting node, but we can we nevertheless rewrite this problem in terms of several instances of the original problem, computed on a subgraph of the original graph.

In particular, let S(i, G) be the set of successors of vertex i in graph G -- that is, the set of vertices s such that there is an edge from i to s in G, and let X(G,i) be the subgraph of G formed by removing all edges that start at i. Furthermore, let p_is be the probability associated with edge (i, s) in G.

Pr(C_if(0)|G) = Sum over s in S(i, G) of p_is * Pr(f|s,X(G,i))

In other words, the probability of reaching f from i through G without revisiting i in between is the sum over all successors of i of the product of the probability of reaching s from i in one step with the probability of reaching f from s through G without traversing any edges outbound from i. That applies for all f in G, including i.

Now observe that S(i, G) and all the p_is are knowns, and that the problem of computing Pr(f|s,X(G,i)) is a new, strictly smaller instance of the original problem. Thus, this computation can be performed recursively, and such a recursion is guaranteed to terminate. It may nevertheless take a long time if your graph is complex, and it looks like a naive implementation of this recursive approach would scale exponentially in the number of nodes. There are ways you could speed the computation in exchange for higher memory usage (i.e. memoization).

---

There are likely other possibilities as well. For example, I'm suspicious that there may be a bottom-up dynamic programming approach to a solution, but I haven't been able to convince myself that loops in the graph don't present an insurmountable problem there.

Answer 2

Problem Clarification

The input data is a set of m rows of n columns of probabilities, essentially an m by n matrix, where m = n = number of vertices on a directed graph. Rows are edge origins and columns are edge destinations. We will, on the bases of the mention of cycles in the question, that the graph is cyclic, that at least one cycle exists in the graph.

Let's define the starting vertex as s. Let's also define a terminal vertex as a vertex for which there are no exiting edges and the set of them as set T with size z. Therefore we have z sets of routes from s to a vertex in T, and the set sizes may be infinite due to cycles ¹. In such a scenario, one cannot conclude that a terminal vertex will be reached in an arbitrarily large number of steps.

In the input data, probabilities for rows that correspond with vertices not in T are normalized to total to 1.0. We shall assume the Markov property, that the probabilities at each vertex do not vary with time. This precludes the use of probability to prioritize routes in a graph search ².

Finite math texts sometimes name example problems similar to this question as Drunken Random Walks to underscore the fact that the walker forgets the past, referring to the memory-free nature of Markovian chains.

Applying Probability to Routes

The probability of arriving at a terminal vertex can be expressed as an infinite series sum of products.

P_t = lim _{s -> ∞} Σ ∏ P_{i, j},

where s is the step index, t is a terminal vertex index, i ∈ [1 .. m] and j ∈ [1 .. n]

Reduction

When two or more cycles intersect (sharing one or more vertices), analysis is complicated by an infinite set of patterns involving them. It appears, after some analysis and review of relevant academic work, that arriving at an accurate set of terminal vertex arrival probabilities with today's mathematical tools may best be accomplished with a converging algorithm.

A few initial reductions are possible.

1. The first consideration is to enumerate the destination vertex, which is easy since the corresponding rows have probabilities of zero.

2. The next consideration is to differentiate any further reductions from what the academic literature calls irreducible sub-graphs. The below depth first algorithm remembers which vertices have already been visited while constructing a potential route, so it can be easily retrofitted to identify which vertices are involved in cycles. However it is recommended to use existing well tested, peer reviewed graph libraries to identify and characterize sub-graphs as irreducible.

Mathematical reduction of irreducible portions of the graph may or may not be plausible. Consider starting vertex A and sole terminating vertex B in the graph represented as {A->C, C->A, A->D, D->A, C->D, D->C, C->B, D->B}.

Although one can reduce the graph to probability relations absent of cycles through vertex A, the vertex A cannot be removed for further reduction without either modifying probabilities of vertices exiting C and D or allowing both totals of probabilities of edges exiting C and D to be less than 1.0.

Convergent Breadth First Traversal

A breadth first traversal that ignores revisiting and allows cycles can iterate step index s, not to some fixed s_max but to some sufficiently stable and accurate point in a convergent trend. This approach is especially called for if cycles overlap creating bifurcations in the simpler periodicity caused by a single cycle.

Σ P_s^{Δ s}.

For the establishment of a reasonable convergence as s increases, one must determine the desired accuracy as a criteria for completing convergence algorithm and a metric for measuring accuracy by looking at longer term trends in results at all terminal vertices. It may be important to provide a criteria where the sum of terminal vertex probabilities is close to unity in conjunction with the trend convergence metric, as both a sanity check and an accuracy criteria. Practically, four convergence criteria may be necessary ³.

1. Per terminal vertex probability trend convergence delta
2. Average probability trend convergence delta
3. Convergence of total probability on unity
4. Total number of steps (to cap depth for practical computing reasons)

Even beyond these four, the program may need to contain a trap for an interrupt that permits the writing and subsequent examination of output after a long wait without the satisfying of all four above criteria.

An Example Cycle Resistant Depth First Algorithm

There are more efficient algorithms than the following one, but it is fairly comprehensible, it compiles without warning with C++ -Wall, and it produces the desired output for all finite and legitimate directed graphs and start and destination vertices possible ⁴. It is easy to load a matrix in the form given in the question using the addEdge method ⁵.

#include <iostream>
#include <list>

class DirectedGraph {

    private:
        int miNodes;
        std::list<int> * mnpEdges;
        bool * mpVisitedFlags;

    private:
        void initAlreadyVisited() {
            for (int i = 0; i < miNodes; ++ i)
                mpVisitedFlags[i] = false;
        }

        void recurse(int iCurrent, int iDestination,
               int route[], int index,
               std::list<std::list<int> *> * pnai) {

            mpVisitedFlags[iCurrent] = true;
            route[index ++] = iCurrent;

            if (iCurrent == iDestination) {
                auto pni = new std::list<int>;
                for (int i = 0; i < index; ++ i)
                    pni->push_back(route[i]);
                pnai->push_back(pni);

            } else {
                auto it = mnpEdges[iCurrent].begin();
                auto itBeyond = mnpEdges[iCurrent].end();
                while (it != itBeyond) {
                    if (! mpVisitedFlags[* it])
                        recurse(* it, iDestination,
                                route, index, pnai);
                    ++ it;
                }
            }

            -- index;
            mpVisitedFlags[iCurrent] = false;
        } 

    public:
        DirectedGraph(int iNodes) {
            miNodes = iNodes;
            mnpEdges = new std::list<int>[iNodes];
            mpVisitedFlags = new bool[iNodes];
        }

        ~DirectedGraph() {
            delete mpVisitedFlags;
        }

        void addEdge(int u, int v) {
            mnpEdges[u].push_back(v);
        }

        std::list<std::list<int> *> * findRoutes(int iStart,
                int iDestination) {
            initAlreadyVisited();
            auto route = new int[miNodes];
            auto pnpi = new std::list<std::list<int> *>();
            recurse(iStart, iDestination, route, 0, pnpi);
            delete route;
            return pnpi;
        }
};

int main() {

    DirectedGraph dg(5);

    dg.addEdge(0, 1);
    dg.addEdge(0, 2);
    dg.addEdge(0, 3);
    dg.addEdge(1, 3);
    dg.addEdge(1, 4);
    dg.addEdge(2, 0);
    dg.addEdge(2, 1);
    dg.addEdge(4, 1);
    dg.addEdge(4, 3);

    int startingNode = 2;
    int destinationNode = 3;

    auto pnai = dg.findRoutes(startingNode, destinationNode);

    std::cout
            << "Unique routes from "
            << startingNode
            << " to "
            << destinationNode
            << std::endl
            << std::endl;

    bool bFirst;
    std::list<int> * pi;
    auto it = pnai->begin();
    auto itBeyond = pnai->end();
    std::list<int>::iterator itInner;
    std::list<int>::iterator itInnerBeyond;
    while (it != itBeyond) {
        bFirst = true;
        pi = * it ++;
        itInner = pi->begin();
        itInnerBeyond = pi->end();
        while (itInner != itInnerBeyond) {
            if (bFirst)
                bFirst = false;
            else
                std::cout << ' ';
            std::cout << (* itInner ++);
        }
        std::cout << std::endl;
        delete pi;
    }

    delete pnai;

    return 0;
}

---

Notes

[1] Improperly handled cycles in a directed graph algorithm will hang in an infinite loop. (Note the trivial case where the number of routes from A to B for the directed graph represented as {A->B, B->A} is infinity.)

[2] Probabilities are sometimes used to reduce the CPU cycle cost of a search. Probabilities, in that strategy, are input values for meta rules in a priority queue to reduce the computational challenge very tedious searches (even for a computer). The early literature in production systems termed the exponential character of unguided large searches Combinatory Explosions.

[3] It may be practically necessary to detect breadth first probability trend at each vertex and specify satisfactory convergence in terms of four criteria

1. Δ(Σ∏P)_t <= Δ_max ∀ t
2. Σ_t=0^T Δ(Σ∏P)_t / T <= Δ_ave
3. |Σ Σ∏P - 1| <= u_max, where u is the maximum allowable deviation from unity for the sum of final probabilities
4. s < S_max

[4] Provided there are enough computing resources available to support the data structures and ample time to arrive at an answer for the given computing system speed.

[5] You can load DirectedGraph dg(7) with the input data using two loops nested to iterate through the rows and columns enumerated in the question. The body of the inner loop would simply be a conditional edge addition.

if (prob != 0) dg.addEdge(i, j);

Variable prob is P _m,n. Route existence is only concerned with zero/nonzero status.