The calculation of modularity seems different with what I've learned

[JamesWZM]

Here is the modularity I know.

And here is your code:

def modularity_matrix(adj_matrix : np.ndarray) -> np.ndarray:
    k_i = np.expand_dims(adj_matrix.sum(axis=1), axis=1)
    k_j = k_i.T
    norm = 1 / k_i.sum()
    K = norm * np.matmul(k_i, k_j)

    return norm * (adj_matrix - K)

I think you want "K" as 1/2m * kikj，and "norm" is the 1/2m，but where is the 1/2?

[JamesWZM]

Ah..... sorry for the wrong expression figure, this one here is the correct one.

[JamesWZM]

I'm sorry for my mistake. sum(k_i) is the 2m in the expression. But there's another problem
I've tried many modularity calculation functions and their answer is different from yours, which made me believe that your function has some problem.
Eventually, I found it, not in the "modularity_matrix" function, but here:

def modularity(mod_matrix : np.ndarray, communities : list) -> float:
    C = np.zeros_like(mod_matrix)
    for community in communities:
        for i, j in combinations(community, 2):
            C[i, j] = 1.0
            C[j, i] = 1.0

    return np.tril(np.multiply(mod_matrix, C), 0).sum()

The problem is the diagonal elements of "C" here is still 0, and your input array is modularity_matrix, not adj_matrix.
So you discarded part of your data here.
And also, you can't use np.tril here.

This function should be:

def modularity(mod_matrix : np.ndarray, communities : list) -> float:
    C = np.zeros_like(mod_matrix)
    for community in communities:
        for i, j in combinations(community, 2):
            C[i, j] = 1.0
            C[j, i] = 1.0
    C[np.diag_indices_from(C)] = 1.0
    return np.multiply(mod_matrix, C).sum()