Compute the Betti numbers of a graph

betti.py

def betti(G, C=None, verbose=False):
    # G is a networkx graph
    # C is networkx.find_cliques(G)

    # RA, 2017-11-03, CC-BY-4.0

    # Ref: 
    # A. Zomorodian, Computational topology (Notes), 2009
    # http://www.ams.org/meetings/short-courses/zomorodian-notes.pdf

    import itertools
    import numpy as np
    import networkx as nx
    from scipy.sparse import lil_matrix

    def DIAGNOSTIC(*params):
        if verbose: print(*params)

    #
    # 1. Prepare maximal cliques
    #

    # If not provided, compute maximal cliques
    if (C is None): C = nx.find_cliques(G)

    # Sort each clique, make sure it's a tuple
    C = [tuple(sorted(c)) for c in C]

    DIAGNOSTIC("Number of maximal cliques: {} ({}M)".format(len(C), round(len(C) / 1e6)))

    #
    # 2. Enumerate all simplices
    #

    # S[k] will hold all k-simplices
    # S[k][s] is the ID of simplex s
    S = []
    for k in range(0, max(len(s) for s in C)):
        # Get all (k+1)-cliques, i.e. k-simplices, from max cliques mc
        Sk = set(c for mc in C for c in itertools.combinations(mc, k + 1))
        # Check that each simplex is in increasing order
        assert (all((list(s) == sorted(s)) for s in Sk))
        # Assign an ID to each simplex, in lexicographic order
        S.append(dict(zip(sorted(Sk), range(0, len(Sk)))))

    for (k, Sk) in enumerate(S):
        DIAGNOSTIC("Number of {}-simplices: {}".format(k, len(Sk)))

    # Euler characteristic
    ec = sum(((-1) ** k * len(S[k])) for k in range(0, len(S)))

    DIAGNOSTIC("Euler characteristic:", ec)

    #
    # 3. Construct the boundary operator
    #

    # D[k] is the boundary operator 
    #      from the k complex 
    #      to the k-1 complex
    D = [None for _ in S]

    # D[0] is the zero matrix
    D[0] = lil_matrix((1, G.number_of_nodes()))

    # Construct D[1], D[2], ...
    for k in range(1, len(S)):
        D[k] = lil_matrix((len(S[k - 1]), len(S[k])))
        SIGN = np.asmatrix([(-1) ** i for i in range(0, k + 1)]).transpose()

        for (ks, j) in S[k].items():
            # Indices of all (k-1)-subsimplices s of the k-simplex ks
            I = [S[k - 1][s] for s in sorted(itertools.combinations(ks, k))]
            D[k][I, j] = SIGN

    for (k, d) in enumerate(D):
        DIAGNOSTIC("D[{}] has shape {}".format(k, d.shape))

    # Check that D[k-1] * D[k] is zero
    assert (all((0 == np.dot(D[k - 1], D[k]).count_nonzero()) for k in range(1, len(D))))

    #
    # 4. Compute rank and dimker of the boundary operators
    #

    # Rank and dimker
    rk = [np.linalg.matrix_rank(d.todense()) for d in D]
    ns = [(d.shape[1] - rk[n]) for (n, d) in enumerate(D)]

    DIAGNOSTIC("rk:", rk)
    DIAGNOSTIC("ns:", ns)

    #
    # 5. Infer the Betti numbers
    #

    # Betti numbers
    # B[0] is the number of connected components
    B = [(n - r) for (n, r) in zip(ns[:-1], rk[1:])]

    DIAGNOSTIC("B:", ns)

    ec_alt = sum(((-1) ** k * B[k]) for k in range(0, len(B)))
    DIAGNOSTIC("Euler characteristic (from Betti numbers):", ec_alt)

    # Check: Euler-Poincare formula
    assert (ec == ec_alt)

    return B


def test1():
    import networkx as nx
    G = nx.Graph()
    G.add_edges_from([(0, 1), (1, 2), (0, 2), (0, 3), (1, 3), (2, 3)])
    G.add_edges_from([(4, 5), (5, 6), (6, 7), (7, 4)])

    # Diagnostic
    print("Number of nodes: {}, edges: {}".format(G.number_of_nodes(), G.number_of_edges()))
    print("Betti numbers:", betti(G, verbose=True))


def test2(n):
    import networkx as nx
    G = nx.cycle_graph(n)
    print(list(G.edges))
    print("Betti numbers:", betti(G, verbose=True))


if __name__ == '__main__':
    test1()

    # this is OK
    test2(3)

    # this fails, probably due to wrong
    # computation of the Euler characteristic
    test2(4)

Dear Numpde

How to limit Betti number only dim=0,1,2. Otherwise it might consume too large of memory.Tks.

with best regards,
Peter

Hi Peter. I suppose you need to pass C containing nodes, edges and triangles. Use enumerate_all_cliques for that.

Best wishes ~

Dear Numpde

Thanks. Solved my problem.

Best regards,
Peter