What is the best known transitive closure algorithm for a directed graph?

In terms of runtime, what is the best known transitive closure algorithm for directed graphs?

I am curre开发者_高级运维ntly using Warshall's algorithm but its O(n^3). Although, due to the graph representation my implementation does slightly better (instead of checking all edges, it only checks all out going edges). Is there any transitive closure algorithm which is better than this? In particular, is there anything specifically for shared memory multi-threaded architectures?

This paper discusses the performance of various transitive closure algorithms:

http://www.vldb.org/conf/1988/P382.PDF

One interesting idea from the paper is to avoid recomputing the entire closure as the graph changes.

There is also this page by Esko Nuutila, which lists a couple of more recent algorithms:

http://www.cs.hut.fi/~enu/tc.html

His PhD thesis listed on that page may be the best place to start:

http://www.cs.hut.fi/~enu/thesis.html

From that page:

The experiments also indicate that with the interval representation and the new algorithms, the transitive closure can be computed typically in time linear to the size of the input graph.

The Algorithm Design manual has some useful information. Key points:

• Transitive closure is as difficult as matrix multiplication; so the best known bound is the Coppersmith–Winograd algorithm which runs in O(n^2.376), but in practice it's probably not worthwhile to use matrix multiplication algorithms.
• For a heuristic speedup, calculate strongly connected components first.

Surprisingly I was unable to find any implementations of the STACK_TC algorithm described by Esko Nuutila (linked by AmigoNico in the other answer).

So I wrote my own simple implementation, in C++. It differs slightly from the original algorithm, see comments in the code for explanation.

It successfully passes a few tests I tried, but readers are encouraged to test it more, and to verify it against the original paper. The code is probably underoptimized.

struct TransitiveClosure
{
    // The nodes of the graph are grouped into 'components'.
    // In a component, each node is reachable (possibly indirectly) from every other node.
    // If the number of components is same as the number of nodes, your graph has no cycles.
    // Otherwise there will be less components.
    // The components form a graph called 'condensation graph', which is always acyclic.
    // The components are numbered in a way that 'B is reachable from A' implies `B <= A`.

    struct Node
    {
        // An arbitrarily selected node in the same component. Same for all nodes in this component.
        std::size_t root = -1;
        // Index if the component this node belongs to.
        std::size_t comp = -1;
    };
    std::vector<Node> nodes; // Size is the number of nodes, which was passed in the argument.

    struct Component
    {
        // Nodes that are a part of this component.
        std::vector<std::size_t> nodes;
        // Which components are reachable (possibly indirectly) from this one.
        // Unordered, but has no duplicates. May or may not contain itself.
        std::vector<std::size_t> next;
        // A convenicene array.
        // `next_contains[i]` is 1 if and only if `next` contains `i`.
        // Some trailing zeroes might be missing, check the size before accessing it.
        std::vector<unsigned char/*boolean*/> next_contains;
    };
    std::vector<Component> components; // Size is the number of components.
};

[[nodiscard]] TransitiveClosure ComputeTransitiveClosure(std::size_t n, std::function<bool(std::size_t a, std::size_t b)> have_edge_from_to)
{
    // Implementation of the 'STACK_TC' algorithm, described by Esko Nuutila (1995), in
    // 'Efficient Transitive Closure Computation in Large Digraphs'.

    constexpr std::size_t nil = -1;

    TransitiveClosure ret;
    ret.nodes.resize(n);
    std::vector<std::size_t> vstack, cstack; // Vertex and component stacks.
    vstack.reserve(n);
    cstack.reserve(n);

    auto StackTc = [&](auto &StackTc, std::size_t v)
    {
        if (ret.nodes[v].root != nil)
            return; // We already visited `v`.
        ret.nodes[v].root = v;
        ret.nodes[v].comp = nil;
        vstack.push_back(v);
        std::size_t saved_height = cstack.size();
        bool self_loop = false;
        for (std::size_t w = 0; w < n; w++)
        {
            if (!have_edge_from_to(v, w))
                continue;
            if (v == w)
            {
                self_loop = true;
            }
            else
            {
                StackTc(StackTc, w);
                if (ret.nodes[w].comp == nil)
                    ret.nodes[v].root = std::min(ret.nodes[v].root, ret.nodes[w].root);
                else
                    cstack.push_back(ret.nodes[w].comp);

                // The paper that this is based on had an extra condition on this last `else` branch,
                // which I wasn't able to understand: `if (v,w) is not a forward edge`.
                // However! Ivo Gabe de Wolff (2019) in "Higher ranked region inference for compile-time garbage collection"
                // says that it doesn't affect correctness:
                // > In the loop over the successors, the original algorithm Stack_TC checks whether
                // > an edge is a so called forward edge. We do not perform this check, which may cause
                // > that a component is pushed multiple times to cstack. As duplicates are removed in the
                // > topological sort, these will be removed later on and not cause problems with correctness.
            }
        }
        if (ret.nodes[v].root == v)
        {
            std::size_t c = ret.components.size();
            ret.components.emplace_back();
            TransitiveClosure::Component &this_comp = ret.components.back();

            this_comp.next_contains.assign(ret.components.size(), false); // Sic.

            if (vstack.back() != v || self_loop)
            {
                this_comp.next.push_back(c);
                this_comp.next_contains[c] = true;
            }

            // Topologically sort a part of the component stack.
            std::sort(cstack.begin() + saved_height, cstack.end(), [&comp = ret.components](std::size_t a, std::size_t b) -> bool
            {
                if (b >= comp[a].next_contains.size())
                    return false;
                return comp[a].next_contains[b];
            });
            // Remove duplicates.
            cstack.erase(std::unique(cstack.begin() + saved_height, cstack.end()), cstack.end());

            while (cstack.size() != saved_height)
            {
                std::size_t x = cstack.back();
                cstack.pop_back();
                if (!this_comp.next_contains[x])
                {
                    if (!this_comp.next_contains[x])
                    {
                        this_comp.next.push_back(x);
                        this_comp.next_contains[x] = true;
                    }

                    this_comp.next.reserve(this_comp.next.size() + ret.components[x].next.size());
                    for (std::size_t c : ret.components[x].next)
                    {
                        if (!this_comp.next_contains[c])
                        {
                            this_comp.next.push_back(c);
                            this_comp.next_contains[c] = true;
                        }
                    }
                }
            }

            std::size_t w;
            do
            {
                w = vstack.back();
                vstack.pop_back();
                ret.nodes[w].comp = c;
                this_comp.nodes.push_back(w);
            }
            while (w != v);
        }
    };

    for (std::size_t v = 0; v < n; v++)
        StackTc(StackTc, v);

    return ret;
}

My test-cases (from the same paper):

• Input: (adjacency matrix, Y is edge source, X is edge destination)

  {0,1,0,0,0,0,0,0},
  {0,0,1,1,0,0,0,0},
  {1,0,0,1,0,0,0,0},
  {0,0,0,0,1,1,0,0},
  {0,0,0,0,0,0,1,0},
  {0,0,0,0,1,0,0,1},
  {0,0,0,0,1,0,0,0},
  {0,0,0,0,0,1,0,0},
  

  Output:

  {nodes=[(0,3),(0,3),(0,3),(3,2),(4,0),(5,1),(4,0),(5,1)],components=[{nodes=[6,4],next=[0],next_contains=[1]},{nodes=[7,5],next=[1,0],next_contains=[1,1]},{nodes=[3],next=[0,1],next_contains=[1,1,0]},{nodes=[2,1,0],next=[3,2,0,1],next_contains=[1,1,1,1]}]}
  

• Input:

  //      a b c d e f g h i j
  /* a */{0,1,0,0,0,1,0,1,0,0},
  /* b */{1,0,1,0,0,0,0,0,0,0},
  /* c */{0,1,0,1,0,0,0,0,0,0},
  /* d */{0,0,0,0,1,0,0,0,0,0},
  /* e */{0,0,0,1,0,0,0,0,0,0},
  /* f */{0,0,0,0,0,0,1,0,0,0},
  /* g */{0,0,0,1,0,1,0,0,0,0},
  /* h */{0,0,0,0,0,0,0,0,1,0},
  /* i */{0,0,1,0,1,0,0,1,0,1},
  /* j */{0,0,0,0,0,0,0,0,0,0},
  

  Output:

  {nodes=[(0,3),(0,3),(0,3),(3,0),(3,0),(5,1),(5,1),(0,3),(0,3),(9,2)],components=[{nodes=[4,3],next=[0],next_contains=[1]},{nodes=[6,5],next=[1,0],next_contains=[1,1]},{nodes=[9],next=[],next_contains=[0,0,0]},{nodes=[8,7,2,1,0],next=[3,2,0,1],next_contains=[1,1,1,1]}]}