Test the hypothesis that the expected number of edges of a random connected graph is \$ O(n \log n) \$

Question

Motivation

The most common model for a random graph is the Erdős–Rényi model. However, it does not guarantee the connectedness of the graph. Instead, let's consider the following algorithm (in python-style pseudocode) for generating a random connected graph with \$n\$ nodes:

g = empty graph
g.add_nodes_from(range(n))

while not g.is_connected:
    i, j = random combination of two (distinct) nodes in range(n)
    if {i, j} not in g.edges:
        g.add_edge(i, j)

return g

The graph generated this way is guaranteed to be connected. Now, my intuition tells me that its expected number of edges is of the order \$ O(n \log n) \$, and I want to test my hypothesis in Python. I don't intend to do a rigorous mathematical proof or a comprehensive statistical inference, just some basic graph plotting.

The Codes

In order to know whether a graph is connected, we need a partition structure (i.e. union-find). I first wrote a Partition class in the module partition.py. It uses path compression and union by weights:

# partition.py

class Partition:
    """Implement a partition of a set of items to disjoint subsets (groups) as
    a forest of trees, in which each tree represents a separate group.
    Two trees represent the same group if and only if they have the same root.
    Support union operation of two groups.
    """

    def __init__(self, items):
        items = list(items)

        # parents of every node in the forest
        self._parents = {item: item for item in items}

        # the sizes of the subtrees
        self._weights = {item: 1 for item in items}

    def __len__(self):
        return len(self._parents)

    def __contains__(self, item):
        return item in self._parents

    def __iter__(self):
        yield from self._parents

    def find(self, item):
        """Return the root of the group containing the given item.
        Also reset the parents of all nodes along the path to the root.
        """
        if self._parents[item] == item:
            return item
        else:
            # find the root and recursively set all parents to it
            root = self.find(self._parents[item])
            self._parents[item] = root
            return root

    def union(self, item1, item2):
        """Merge the two groups (if they are disjoint) containing
        the two given items.
        """
        root1 = self.find(item1)
        root2 = self.find(item2)

        if root1 != root2:
            if self._weights[root1] < self._weights[root2]:
                # swap two roots so that root1 becomes heavier
                root1, root2 = root2, root1

            # root1 is heavier, reset parent of root2 to root1
            # also update the weight of the tree at root1
            self._parents[root2] = root1
            self._weights[root1] += self._weights[root2]

    @property
    def is_single_group(self):
        """Return true if all items are contained in a single group."""
        # we just need one item, any item is ok
        item = next(iter(self))

        # group size is the weight of the root
        group_size = self._weights[self.find(item)]
        return group_size == len(self)

Next, since we are only interested in the number of edges, we don't actually need to explicitly construct any graph object. The following function implicitly generates a random connected graph and return its number of edges:

import random
from partition import Partition

def connected_edge_count(n):
    """Implicitly generate a random connected graph and return its number of edges."""
    edges = set()
    forest = Partition(range(n))

    # each time we join two nodes we merge the two groups containing them
    # the graph is connected iff the forest of nodes form a single group
    while not forest.is_single_group:
        start = random.randrange(n)
        end = random.randrange(n)

        # we don't bother to check whether the edge already exists
        if start != end:
            forest.union(start, end)
            edge = frozenset({start, end})
            edges.add(edge)

    return len(edges)

We then estimate the expected number of edges for a given \$n\$:

def mean_edge_count(n, sample_size):
    """Compute the sample mean of numbers of edges in a sample of
    random connected graphs with n nodes.
    """
    total = sum(connected_edge_count(n) for _ in range(sample_size))
    return total / sample_size

Now, we can plot the expected numbers of edges against \$ n \log n \$ for different values of \$n\$:

from math import log
import matplotlib.pyplot as plt

def plt_mean_vs_nlogn(nlist, sample_size):
    """Plot the expected numbers of edges against n * log(n) for
    a given list of values of n, where n is the number of nodes.
    """
    x_values = [n * log(n) for n in nlist]
    y_values = [mean_edge_count(n, sample_size) for n in nlist]
    plt.plot(x_values, y_values, '.')

Finally, when we called plt_mean_vs_nlogn(range(10, 1001, 10), sample_size=100), we got:

The plot seems very close to a straight line, supporting my hypothesis.

Questions and ideas for future work

1. My program is slow! It took me 90 seconds to run plt_mean_vs_nlogn(range(10, 1001, 10), sample_size=100). How can I improve the performance?
2. What other improvement can I make on my codes?
3. An idea for future work: do a linear regression on the data. A high coefficient of determination would support my hypothesis. Also find out the coefficient of \$ n \log n \$.
4. Any other idea for testing my hypothesis programmatically?

Answer 1

"My program is slow!"

You want an estimate for \$P\$ different graph-sizes, each of which is the average of \$S\$ samples of connected_edge_count. We assume connected_edge_count will run through it's while loop \$n\log n\$ times (approximately). What's the asymptotic complexity of Partition.find()? I'll wildly guess it's \$\log n\$. So taking \$N\$ as the maximum requested \$n\$, your overall program is \$O(P S N (\log N)^2)\$.

So broadly speaking, there's just a lot of work to do. Local improvements to the implementation details can help, but I think your biggest problem (at least until you start increasing \$n\$) is \$S\$. 100 is way too big. Playing around with some values, 15 seems to give somewhat stable results, although possibly you'll want more as you deal with larger \$n\$ values.

On the other hand, how often are you planning on running this? Ninety seconds isn't really that long. It feels like a long time when you're trying to iterate on the program. So one thing you might want to work on is the way the functions are nested. Rather than having each function in the stack call the next, let them take the prior result as an argument. This way you'll have better access to intermediate results, and won't have to re-run everything every time.

I spent some time squishing around parts of the code to make sure I understood it, and then because I couldn't get the details out of my head. I haven't checked if it's faster or now, mostly it's just denser. For an academic POC, it goes up to 10K without any problems. (My main takes about three minutes to run. I still can't get connected_edge_count(10 * 1000 * 1000) to work; it crashes after a few minutes.) I'll post my version below in case there are any differences in it you find useful.

"What other improvement can I make on my codes?"
All the usual stuff. Better variable names, less mutation of state and variables, type-hints. Once I got a sense of what your code did I quite liked it; the tree system is clever. (But is it accurate? How do you know? If you're hoping to publish results, adding a few unit tests isn't going to be good enough.)

In your comments you talked about not needing to build an explicit graph; you claimed to do it virtually. But notice that you do need to keep track of all the edges so that you can count them.

Because performance is an issue, and we want to be able to handle large numbers of items, I made some optimizations that may make the code harder to read. For example, for the task at hand a List[int] (array) can serve the purpose of a Dict[int, int] with a lot less machine overhead. But it ties you representing your nodes as ints.

As for further research steps, it depends on your goals. My intuition is that this kind of sampling may be an easy way of checking if your hypothesis is viable, and you've done that. If you want to prove it, then you need to actually prove it. Maybe a programmatic proof system like agda or coq can help, maybe not; I haven't learned them yet!

import matplotlib
matplotlib.use('TkAgg')

from itertools import count, dropwhile, takewhile
import random
from math import exp, log
import matplotlib.pyplot as plt
from scipy.special import lambertw
from typing import Callable, Dict, Iterable, List, NamedTuple, Set, Tuple

from time import sleep


class Partition(NamedTuple):
    parents: List[int]
    weights: List[int]
    edges: Set[Tuple[int, int]]  # The tuple members must be storted! frozensets would be harder to screw up, and maybe slightly faster, but they take more ram, which I think is the bottleneck.

    @staticmethod
    def atomic(node_count: int):
        return Partition(
            parents=list(range(node_count)),
            weights=[1 for node in range(node_count)],
            edges=set()
        )

    def _node_to_str(self, node: int) -> str:
        if not node < len(self.parents):
            raise Exception(f"{node} is not in the range 0 - {len(self.parents)}.")
        return "{n}: <{c}>".format(
            n=node,
            c=", ".join(self._node_to_str(n) for (n, p) in enumerate(self.parents) if p == node and n != node)
        )

    def display(self) -> str:
        if 100 < len(self.parents):
            raise NotImplementedError("Refusing to pretty-print a large forest.")
        return "\n".join(self._node_to_str(n) for (n, p) in enumerate(self.parents) if p == n)

    def find_root(self, item: int) -> int:
        parent = self.parents[item]
        if parent == item:
            return item
        else:  # find the root and recursively set all parents to it
            root = self.find_root(parent)
            self.parents[item] = root
            return root

    def add_edge(self, item1: int, item2: int) -> int:
        """returns the number of edges added to the graph (1, or 0 if the edge was already there)"""
        edge = (item1, item2) if item1 < item2 else (item2, item1)
        if edge in self.edges:
            return 0
        else:
            self.edges.add(edge)
            root1 = self.find_root(item1)
            root2 = self.find_root(item2)
            if root1 != root2:
                weight1 = self.weights[root1]
                weight2 = self.weights[root2]
                heavier, lighter, lesser_weight = (root2, root1, weight1) if weight1 < weight2 else (root1, root2, weight2)
                self.parents[lighter] = heavier  # reset parent of lighter to heavier
                self.weights[heavier] += lesser_weight  # also update the weight of the tree the heavier node
            return 1

    def is_single_group(self) -> bool:
        # we can start with any node for the task at hand
        return len(self.parents) == self.weights[self.find_root(self.parents[0])]


def connected_edge_count(n: int) -> int:
    forest = Partition.atomic(n)
    nodes = range(n)  # not the _real_ nodes, just an external thing to sample from.
    while not forest.is_single_group():
        edge = random.sample(nodes, 2)
        forest.add_edge(*edge)
    return len(forest.edges)


def mean_of(trial: Callable[..., int], *trial_args, sample_size: int, **trial_kwargs) -> float:
    return sum(trial(*trial_args, **trial_kwargs) for _ in range(sample_size)) / sample_size


def nlogn(x):
    return x * log(x)


def inverse(x):
    return abs(x / lambertw(x))


def plt_vs_nlogn(*samples: Tuple[int, float]):
    x_values = [nlogn(n) for (n, v) in samples]
    plt.xlabel("n⋅log(n)")
    y_values = [v for (n, v) in samples]
    plt.ylabel("mean edges to connect n-graph")
    plt.plot(x_values, y_values, '.')


def nlogn_range(start: int, stop: int, starting_step: float = 100) -> Iterable[int]:
    """This is rediculious overkill."""
    return [
        int(inverse(x))
        for x
        in takewhile(lambda _x: inverse(_x) < stop,
                     dropwhile(lambda _x: inverse(_x) < start,
                               count(1, nlogn(starting_step))))
    ]


def main():
    ns = list(nlogn_range(10, 10 * 1000, 500))
    results = {
        n: mean_of(
            connected_edge_count,
            n,
            sample_size=int(5 * log(n))
        )
        for n in ns
    }
    plt_vs_nlogn(*results.items())


def agrees_with_original(i: int) -> bool:
    import cr_243594_original
    mine = mean_of(connected_edge_count, i, sample_size=i)
    theirs = cr_243594_original.mean_edge_count(i, i)
    print(mine)
    print(theirs)
    return abs(mine - theirs) < (i/10)  # this is totally made up and prone to failure because of the randomness.

def verbose_random_add(tree: Partition) -> None:
    edge = random.sample(range(len(tree.parents)), 2)
    print(sorted(edge))
    print(sorted(tree.edges))
    print(tree.add_edge(*edge))
    print(tree.display())
    print(tree.is_single_group())
    print(sorted(tree.edges))
    assert all(((a,b) not in tree.edges) or (tree.find_root(a) == tree.find_root(b))
               for a in range(len(tree.parents))
               for b in range(a+1, len(tree.parents)))
    print("===========================================")


assert agrees_with_original(40)


if __name__ == "__main__":
    main()