Counting triangles in a graph by iteratively removing high-degree nodes

Assuming the graph is not directed (ie. G.is_directed() == False), the number of triangles can be efficiently found by finding nodes that are both neighbors of neighbors and direct neighbors of a same node. Pre-computing and pre-filtering the neighbors of nodes so that each triangle is counted only once helps to improve a lot the execution time. Here is the code:

nodeNeighbours = {
    # The filtering of the set ensure each triangle is only computed once
    node: set(n for n in edgeInfos.keys() if n > node)
    for node, edgeInfos in G.adjacency()
}

triangleCount = sum(
    len(neighbours & nodeNeighbours[node2])
    for node1, neighbours in nodeNeighbours.items()
    for node2 in neighbours
)

The above code is about 12 times faster than the original iterative solution on the example graph. And up to 72 times faster on nx.erdos_renyi_graph(15000, 0.005).

A faster solution consist in using Numba and perform the computation with multiple threads. Unfortunately, this method is much more complex (hence this separate answer).

The idea is to first convert the graph into Numpy arrays, then balance the work between the thread, and finally apply the algorithm provided in the alternative answer (using multiple thread and Numpy arrays). The graphs is encoded using four Numpy arrays:

• nodes: contains all the nodes IDs;
• allNeighbours: contains all the neighbours of all nodes;
• neighbourSlices: contains for each node, the slice (start+end index) of the array allNeighbours (so to be able to fetch the neighbours of a node);
• nodeIdToPos: contains the index of the node in nodes based on its ID.

Here is the code:

# Split the nodes so the work (based on the slices) is balanced.
# Returns the (start,stop) slices of the nodes
@nb.njit('(int_[:,::1], int_)')
def splitSlices(neighbourSlices, count):
    n = np.int64(neighbourSlices.shape[0])
    m = neighbourSlices[-1, 1] if neighbourSlices.size > 0 else 0
    workSize = np.empty((count, 2), dtype=np.int_)
    curPos = 0
    for i in range(count):
        for j in range(curPos, n):
            stop = neighbourSlices[j, 1]
            if stop >= m * (i + 1) // count:
                workSize[i, 0] = curPos
                curPos = j + 1
                workSize[i, 1] = curPos
                break
    if count > 0:
        workSize[0, 0] = 0
        workSize[-1, 1] = n
    return workSize


@nb.njit('(int_[::1], int_[:,::1], int_[::1])', parallel=True)
def filterNeighbours(nodes, neighbourSlices, allNeighbours):
    outNeighbourSlices = np.empty(neighbourSlices.shape, dtype=np.int_)
    outAllNeighbours = np.empty(allNeighbours.size//2, dtype=np.int_)

    curPos = 0
    for i in range(nodes.size):
        curNode = nodes[i]
        start, stop = neighbourSlices[i]
        outNeighbourSlices[i, 0] = curPos
        for neighbour in allNeighbours[start:stop]:
            if neighbour > curNode:
                outAllNeighbours[curPos] = neighbour
                curPos += 1
        outNeighbourSlices[i, 1] = curPos

    threadCount = nb.np.ufunc.parallel.get_num_threads()
    nodeSlides = splitSlices(outNeighbourSlices, threadCount)

    for threadId in nb.prange(threadCount):
        start, stop = nodeSlides[threadId]
        for i in range(start, stop):
            start, stop = outNeighbourSlices[i]
            outAllNeighbours[start:stop].sort()

    return (outNeighbourSlices, outAllNeighbours)


@nb.njit('int_(int_[::1], int_[:,::1], int_[::1])', parallel=True)
def computeTriangle(nodes, neighbourSlices, allNeighbours):
    nodeIdToPos = np.empty(np.max(nodes)+1, dtype=np.int_)
    for i, node in enumerate(nodes):
        nodeIdToPos[node] = i

    neighbourSlices, allNeighbours = filterNeighbours(nodes, neighbourSlices, allNeighbours)

    threadCount = nb.np.ufunc.parallel.get_num_threads()
    nodeSlides = splitSlices(neighbourSlices, threadCount)

    s = 0
    for threadId in nb.prange(threadCount):
        start, stop = nodeSlides[threadId]
        for i in range(start, stop):
            node1 = nodes[i]
            start1, stop2 = neighbourSlices[i]
            neighbours1 = allNeighbours[start1:stop2]
            for node2 in neighbours1:
                start2, stop2 = neighbourSlices[nodeIdToPos[node2]]
                neighbours2 = allNeighbours[start2:stop2]
                for node3 in neighbours2:
                    if node3 > node2:
                        foundId = np.searchsorted(neighbours1, node3)
                        s += foundId < neighbours1.size and neighbours1[foundId] == node3
    return s

# Conversion to Numpy arrays

edgeCount = sum(len(neighbours) for _, neighbours in G.adjacency())
nodes = np.fromiter(G.nodes(), dtype=np.int_)
neighbourSlices = np.empty((len(G), 2), dtype=np.int_)
allNeighbours = np.empty(edgeCount, dtype=np.int_)

curPos = 0
for i, nodeInfos in enumerate(G.adjacency()):
    node, edgeInfos = nodeInfos
    neighbours = np.fromiter(edgeInfos.keys(), dtype=np.int_)
    neighbourSlices[i, 0] = curPos
    neighbourSlices[i, 1] = curPos + neighbours.size
    allNeighbours[curPos:curPos+neighbours.size] = neighbours
    curPos += neighbours.size

# Actual computation (very fast)

triangleCount = computeTriangle(nodes, neighbourSlices, allNeighbours)

On my 6-core machine, this solution 45 times faster than the original iterative solution on the example graph and up to 290 times faster on nx.erdos_renyi_graph(15000, 0.005). It is much much faster than other methods when the average degree (N_edges/N_nodes) is big.

Unfortunately, most of the time is spent in the (sequential) conversion of the graph into Numpy arrays. Thus, if you want a faster method (eg. 2~3 times faster), then you certainly need to use a native language implementation (like C or C++). Or to work mainly with Numpy array instead of NetworkX graphs, but this approach is clearly less flexible and complex to archive.