Hopcroft's Algorithm for Minimization of Deterministic Finite State Automaton

Question

Background

Given a deterministic finite state automaton, Hopcroft's algorithm is a $O(n \log(n))$ algorithm for finding an automaton that generates the same language with a minimal number of vertices. We can represent an automaton $\mathcal{A}$ as a directed labeled graph $\mathcal{A}=(Q,F, \Sigma, \delta)$: a set of vertices/states $Q$, a set of fail vertices/states $F\subseteq Q$, an alphabet of edge labels $\Sigma$, and a collection of labeled edges $\delta:Q\times \Sigma\to Q$.

Copying the pseudocode from wikipedia, we generate a partition $P$ of the vertex set $Q$ by repeatedly splitting the current partition, while updating a waiting set $W$ by which we need to split/refine the current partition.

P := {F, Q \ F};
W := {F};
while (W is not empty) do
     choose and remove a set A from W
     for each c in Σ do
          let X be the set of states for which a transition on c leads to a state in A
           for each set Y in P for which X ∩ Y is nonempty and Y \ X is nonempty do
                replace Y in P by the two sets X ∩ Y and Y \ X
                if Y is in W
                    replace Y in W by the same two sets
                else
                    if |X ∩ Y| <= |Y \ X|
                         add X ∩ Y to W
                    else
                         add Y \ X to W
          end;
     end;
end;

I'm having some issues converting the numerous nested loops into clean Mathematica/functorial code. A naive, faithful implementation While[...For[...For[...If[...,....,If[...,...,...]]]]] is just as slow as one would expect.

My implementation

I'm storing graphs as a list of directed edges, e.g.,

g1={{0 -> 0, 0}, {0 -> 0, 1}, {1 -> 2, 0}, {1 -> 3, 1}, 
    {2 -> 2, 0}, {2 -> 3, 1}, {3 -> 0, 1}, {3 -> 1, 0}};

where the vertex $0$ is the sole fail state.

My first helper function creates a table (indexed by the alphabet) of associations that return all vertices which have an edge bearing the appropriate label leading to the desired vertex, i.e., $\delta^{-1}(x, e)$:

incomingVertices[graph_, vertexList_, alphabet_] := 
    Table[
       With[{relevant = Select[graph, #[[2]] == a &]}, 
            Merge[
             {Thread[relevant[[All, 1, 2]] -> relevant[[All, 1, 1]]], 
           Thread[Complement[vertexList, 
          DeleteDuplicates@relevant[[All, 1, 2]]] -> Nothing]},
     Join]
    ],
   {a, alphabet}] 

A bit kludgy, but it only needs to run once.

incomingVertices[g1, VertexList[g1[[All,1]]], Union[g1[[All,2]]]]
(*{<|0 -> {0}, 2 -> {1, 2}, 1 -> {3}, 3 -> {}|>, 
    <|0 -> {0, 3},  3 -> {1, 2}, 1 -> {}, 2 -> {}|>} *)

Edit: Per Daniel Lichtblau and Henrik Schumacher in the comments, it appears that this implementation is $O(n^2)$ or worse

My main function, with the general structure While[...Do[...Table[]; Table[];..]...], is then given by

hopcroft[graph_] := Module[
{alphabet = Union@graph[[All, 2]],
 vertexList = Sort@VertexList@graph[[All, 1]],
 W = {{0}},
 inc, part, partby, newtopartby, inW, notinW, newpartinW, newpartnotinW, newtopartbyinW, newtopartbynotinW, X},
 part = {DeleteCases[vertexList, 0], {0}};
 inc = incomingVertices[graph, vertexList, alphabet];
 Clear[vertexList];
 While[Length@W > 0,
    partby = First@W;
    W = Rest[W];
    Do[X = Union@Flatten[partby /. inc[[index]]];
        inW=Intersection[part, W];
        notinW=Complement[part, inW];
        newpartinW=Table[{Intersection[Y, X], Complement[Y, X]}, {Y, inW}];
        newpartnotinW=Table[{Intersection[Y, X], Complement[Y, X]}, {Y, notinW}];
        part=DeleteCases[Level[newpartinW~Join~newpartnotinW, {-2}], {}];
        newtopartbyinW=Level[W /. Association[Thread[inW -> newpartinW]], {-2}];
        newtopartbynotinW= MinimalBy[#,Length,1]&/@Select[newpartnotinW,FreeQ[#, {}] &];
        W = DeleteCases[Level[newtopartbyinW~Join~newtopartbynotinW, {-2}], {}];,
{index, Length@alphabet}]
   ];
  part]

The general philosophy of my implementation is that we need to handle pieces of the partition differently depending on whether or not they are in the waiting set W, hence the initial inW = ... and notinW=.... Similarly, my prototyping lead to it seeming to be more efficient to just compute all the intersections and complements at once, rather than implementing a bunch of boolean logic to handle each element of partitions independently.

Timing results

For a small graph, this runs quite fast:

RepeatedTiming[hopcroft[g1]]
(*{0.000310, {{0}, {3}, {1, 2}}}*)

However, it starts grinding to a halt for large graphs. For example, the following sort of graph is generally considered a worst case scenario for the algorithm, needing to refine the partition until each vertex is distinct:

bigGraph = With[{n = 10^3, alph = 3}, 
           Flatten[#, 1]&@
           Table[{i -> If[i + j > n, 0, i + j], j}, {i, n}, {j, alph}]];

RepeatedTiming[hopcroft[bigGraph];]
(*{3.663, Null} *)

With n=10^4 instead, I get

Timing[hopcroft[bigGraph];]
(*{362.041, Null}*)

My ideal end usage is even worse, with automata that have $O(10^5)$ vertices and $O(10^6)$ edges.

Questions

1. How can I speed this up? Profiling my code, it seems to spend by far the majority of the time computing an obscene number of intersections. Since this is already such a low-level built in function, I'm not sure what I can do to significantly speed this up.
2. I have access to some clustered computing resources (my desktop has 4 physical cores, free reign on a server with 32 physical cores, scheduled access to a cluster with over 1000 physical cores). However, I don't see a great way to parallelize this algorithm, as the partition will (usually) get updated with each lap of the Do[X = ..., {index, Length@alphabet}] loop. Naively replacing all the Table[..] commands with ParallelTable[...] seems to actually slow down computation, presumably since the individual computations are so fast and we gain a lot of overhead of pushing data to/from the slave cores. It is, however, quite entertaining to watch the ParallelKernelStatus window while this happens! Any thoughts on how we can efficiently distribute this problem?
3. There are implementations of Hopcroft's algorithm implemented in C/C++/Python floating around---in many ways, the algorithm seems built for procedural languages. Worst case, I could always run an external call. However, I have a preference towards staying within Mathematica---especially if we can efficiently parallelize it.

Edits:

Hopcroft's original 1971 paper can be found here, but is unfortunately paywalled. There are several other analyses and descriptions of the algorithm in the literature as well, primarily explanatory or focused on showing that it is indeed $O(n \log(n))$ and $\Omega(n \log(n))$ or analyzing run time comparisons to several other algorithms that have worse theoretical complexity, but can slightly outperform Hopcroft in some cases.

Explanation of My Implementation

I have also tried to present my code with generally the same variable names as in the pseudocode (particularly W, X, and Y). The main change between the pseudocode an my implementation concerns converting the innermost for ... loop with its multiple booleans into just computing all of the intersections / set complements and then updating the partition part and the waiting set W accordingly. Specifically, my code handles these operations with the innermost Do loop:

1. X = Union@Flatten[partby /. inc[[index]]]; Sets X as the set of all vertices that have an edge labeled by alphabet[[index]]. We need to update the current partition part and waiting set W based on the intersections / complements with X of each set Y in part and W
2. inW=Intersection[part, W]; notinW=Complement[part, inW]; breaks the current partition part into two chunks: those that also appear in the waiting list, and those that don't. Rather than using booleans later to handle these separately, I sepearate them out now.
3. newpartinW=Table[{Intersection[Y, X], Complement[Y, X]}, {Y, inW}]; newpartnotinW=Table[{Intersection[Y, X], Complement[Y, X]}, {Y, notinW}]; For each set $Y$ in the partition part, we will eventually replace $Y$ by $Y \cap X$ and $Y \setminus X$. Note that one of these will often be empty the empty list {}, so we will garbage collect out these later.
4. part=DeleteCases[Level[newpartinW~Join~newpartnotinW, {-2}], {}]; returns the new partition. Regardless of whether Y is in the waiting set or W or not, we replace Y in part by $\{ Y \cap X, Y\setminus X\}$, use Level[ ..., {-2}] to flatten out a level, and then remove all the empty lists with DeleteCases[ ..., {}].
5. newtopartbyinW=Level[W /. Association[Thread[inW -> newpartinW]], {-2}]; If a set $Y$ is in the waiting set $W$, then we need to replace $Y$ in $W$ by $Y \cap X$ and $Y \setminus X$. Since inW contains all elements of the partition that are in W and newpartinW contains the image of each under the mapping $Y \mapsto \{ Y \cap X, Y \setminus X \}$, we use Thread[inW -> newpartinW] to get a list of the form {{Y_1 -> {Y_1 \cap X, Y_1 \setminus X}, {Y_2 ->{Y_2 \cap X, Y_2 \setminus X}, ...} and then update W by acting on it by this association. Note that there will again be many empty lists {} (which will later get removed). Also note that if $w \in W$ and $w \not \in$part, then it will be unchanged, i.e., we retain the previous waiting list W and only split the ones that are elements of the current partition.
6. newtopartbynotinW= MinimalBy[#,Length,1]&/@Select[newpartnotinW,FreeQ[#, {}] &]; This line takes care of the bulk of the booleans Else[ If[ ..., ..., Else[...]]]. The Select[newpartnotinW,FreeQ[#, {}] &] picks out only the pairs ${Y_i \cap X, Y_i \setminus X}$ when $Y_i$ is actually split by $X$ (if a set $Y_i$ isn't split by $X$, then we don't need to add anything to W. From this list of non-empty pairs, MinimalBy[#, Length,1]&/@ ... maps over the list of pairs and returns whichever is shorter, with the MinimalBy[...,...,1] ensuring that only one list from each pair is returned in the event of a tie in length.
7. W = DeleteCases[Level[newtopartbyinW~Join~newtopartbynotinW, {-2}], {}]; We redefine the waiting set W to update all the changes we need to make. Note that we will loop back through the Do[..., {index, Length@alphabet}] several more times before stepping out to pick a new set with which to partition by. The new waiting list W contains sublists from three sources:
   • Those that were already in W and are not elements of the partition. These are in newtopartbyinW as they are left unchanged by W/.Associtation[...]
   • Those that were already in W and are elements of the partition. These are in newtopartbyinW as they were updated by the replacement W/.Association[...]
   • Those that were not in W already. These come from newtopartbynotinW. Note that we only needed to add in the shorter of $Y_i \cap X$ and $Y_i \setminus X$, as splitting by a set $X$ is the same as splitting by $X^c$.

We then repeat the Do loop for the next splitting set X, then break out and choose and delete a new set topartby from the waiting set W.

The key observation that I had regarding this level of optimization was that at any stage in the algorithm, the current partition part (or P) has to satisfy containing every element of the vertex list (i.e., Union[part]==VertexList[graph[[All,1]]), and the partition is in fact a partition of the vertex list (i.e., it is composed of non-intersecting lists). As such, updating a set $Y \in P$ has no side-effects on other elements of the partition and updating W based upon the splitting of a set Y has no side-effects on the other elements of W (up to the specific order of the waiting list/ queue W). We can see this most clearly in that the innermost for loop of the pseudocode doesn't require an order with which to split each of the $Y \in P$. Thus, I have done things in the order "split the $Y$'s that are also in W first" and then "split the $Y$'s thare aren't in W" last.

Some Musings Towards Optimization

Here's some ideas I've had, but haven't quite been able to figure out how to connect the dots.

• There does seem to be some research in the literature into the effect of different strategies for processing the waiting list / queue W, i.e., FIFO vs. FILO vs. Shortest first vs. Longest first, etc. While it does seem to make some difference, the effect is a bit underwhelming (<5% gains under optimal conditions IIRC), so I've elected to just take topartby=First[W]; W=Rest[W]; out of convenience.
• If an element of the partition gets split apart enough to the point that it only contains a single vertex (i.e., {..., {v_i}, ....}), then there's no need to keep keep recomputing $\{ \{v_i\} \cap X, \{v_i\} \setminus X \}$, as one of them will trivially be empty. So perhaps some logic to set such sets aside and treat them differently? Note that singleton sets in the waiting list W are fine and necessary.
• The current structure takes a set X = Union@Flatten@(partby/.inc[[index]]), updates the partition, updates the waiting list, and then loops back around to take X=Union@Flatten(partby/.inc[[index + 1]]. This seems like its perhaps the prone to effective parallelization, sending each core off to compute the new partition and waiting list for its given index. I.e., something like ParallelTable[...., {index, Length@alphabet}] instead of Do[...]. I haven't been able to envision how to combine back these results without computing a whole mess of pairwise intersections/complements. Similarly, the idea of distributing an element of the waiting list to each core feels productive, until its time to combine all of the results back into a single partition and single waiting list.
• All the constant reassignments W=... and part=... are generally slow, as opposed to some faster data structure that doesn't require copying the entire list each time.
• While Intersection and Complement are already fast built-ins, perhaps there's savings to be had by using a compiled function which knows that it will only ever act integer lists? Similarly, $Y \setminus X = Y \setminus (Y \cap X)$. Since $|Y \cap X| \leq | Y|$ and presumably Complement[Y, #] will be faster with a shorter list in the second component, I've tried using functions along the lines of With[{YcapX = Intersection[Y,X]},{YcapX, Complement[Y,YcapX]}]. Unfortunately, this seemed to actually slow things down, presumably due to needing to store that intermediate value in memory.