Newman algorithm yielding different result to what is given in his paper

Question

Summary

I am trying to implement Newman's algorithm for community detection, outlined in this paper. I am testing my implementation against one of the datasets used in that paper to benchmark the algorithm, and I am getting slighty different, and less optimal results.

2 nodes are being placed in the wrong group, which reduces the modularity. I can pinpoint the exact spot at which it "goes wrong" (i've marked a place to put a breakpoint in the code), but i'm not sure how to fix it, or why my approach is wrong. Code below.

Further detail

What's going wrong

My implementation of the algorithm places nodes $1$ and $12$ into the wrong community. The Zachary Karate network is well studied, and other sources (which achieve a similar modularity score as Newman in the paper) have the clustering depicted below. My clustering is added for contrast.

What i've tried

Below is my Python implementation of the algorithm. I originally did this in MATLAB (also running in Octave), both gave the exact same result, which makes me think i've made some fundamental flaw that i'm just unable to see right now.

I also tried using variable precision arithmetic in MATLAB in case this was a round off error, to no avail.

Python code

The implementation goes noticeably wrong when trying to maximise the modularity of the nodes [ 1 2 3 4 5 6 7 8 11 12 13 14 17 18 20 22]. I've purposely left the check for this group in the code so it is easy to set a break point on the correct recursion step.

At this step, the eigenvector components corresponding to nodes $1$ and $12$ (indices 0 and 9 respectively) should have the same sign as those corresponding to nodes $2,3,4,8,13,14,18,20,22$ in leading_eigen_vector, but they do not, instead they have the same signs as $5,6,7,11,17$, which causes them to be placed in the wrong group.

Note I don't calculate the exact modularity score here, but in my MATLAB implementation, this gives a modularity of $Q = 0.3934$, whereas Newman achieves a modularity for this network of $Q = 0.419$. I have also tried using the $\Delta Q$ from the paper to determine when a division is good or not, and the issue persists. If I manually move the 2 incorrectly placed nodes, I achieve the same modularity as Newman.

import numpy as np
from numpy import linalg as LA
import sys
np.set_printoptions(threshold=sys.maxsize)

def communities(B, category, globals):
    print(globals + 1) # debugging code - globals are the nodes we are looking at for this step 

    I = np.identity(B.shape[0])
    B_transpose = B.transpose()

    # kronecker_sum calculates the kronecker delta * sum of B rows (from equation 6)
    kronecker_sum = np.multiply( I , 
                          np.sum(B_transpose, axis = 1).reshape(B.shape[0],1) # sum up the transpose of B, and turn it into a column vector for the next step
                        )

    # Compute equation 6
    Bg = np.subtract(B, kronecker_sum)

    eigen_values, eigen_vectors = LA.eig(Bg)

    # Find the most positive eigenvalue, and the corresponding eigenvector
    leading_eigen_value = np.amax(eigen_values)
    index_of_lead = np.where(eigen_values == leading_eigen_value)
    leading_eigen_vector = eigen_vectors[:, [index_of_lead]] # extract the column vector representing the leading eigenvector

    if np.array_equal(globals + 1, np.array([1, 2, 3, 4, 5, 6, 7, 8, 11, 12, 13, 14, 17, 18, 20, 22])):
        # indices 0 and 9 of leading_eigen_vector will be negative, they should be positive to place nodes 1 and 12 into the correct group
        # that would maximise modularity
        place_break_point_here = True

    # membership vector (place network nodes in 1 group if the same eignevector index is geq to 0, else put into a different group)
    s = np.where(leading_eigen_vector >= 0, 1, -1)

    if (leading_eigen_value < 0.1):
        labels = np.full((1, B.shape[0]), category)
        category = category + 1
        return [labels, category]
    else:
        # node indices in Bg that correspond to the first group
        left_indices  = np.array([elem[0] for elem in np.argwhere(s ==  1)])

        # node indices in Bg that correspond to the second group
        right_indices = np.array([elem[0] for elem in np.argwhere(s == -1)])

        # Elements in B corresponding to nodes in our first and second groups respectively
        left_B =  B[np.ix_(left_indices,left_indices)]   
        right_B = B[np.ix_(right_indices,right_indices)]

        # recurse on our group, try and split them up further
        [left, category]  = communities(left_B, category, globals[np.ix_(left_indices)])   
        [right, category] = communities(right_B, category, globals[np.ix_(right_indices)])

        labelled_vertices = np.zeros(max(left.shape) + max(right.shape)) # allocate an array of the correct size to put our labelled nodes in
        labelled_vertices[np.ix_(left_indices)] = left
        labelled_vertices[np.ix_(right_indices)] = right      
        return [labelled_vertices, category]

# Adjacency matrix from Zachary's Karate dataset
A = np.array([
    [0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0],
    [1, 0, 1, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0],
    [1, 1, 0, 1, 0, 0, 0, 1, 1, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 1, 0],
    [1, 1, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
    [1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
    [1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
    [1, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
    [1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
    [1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1],
    [0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1],
    [1, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
    [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
    [1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
    [1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1],
    [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1],
    [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1],
    [0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
    [1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
    [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1],
    [1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1],
    [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1],
    [1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
    [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1],
    [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 1, 1],
    [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 1, 0, 0],
    [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0],
    [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1],
    [0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1],
    [0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1],
    [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1],
    [0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1],
    [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 1, 1],
    [0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 0, 0, 1, 0, 1, 0, 1, 1, 0, 0, 0, 0, 0, 1, 1, 1, 0, 1],
    [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 1, 1, 1, 0, 1, 1, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0], 
])

degrees = np.sum(A, axis = 1).reshape(A.shape[0],1)
m = np.sum(degrees)/2
K = np.outer(degrees, degrees.transpose()[0])
B = np.subtract(A, K/(2*m))
[labelled_vertices, label] = communities(B, 0, np.arange(A.shape[0]))

A plain English description of the Algorithm

At each pass of maximising the modularity of the network, we compute equation 6 in the paper. This will yield a matrix, for which we compute the eigenvectors and eigenvalues. We look at the eigenvector corresponding with the most positive eigenvalue. We look at each entry of this eigenvector: if the entry is greater than or equal to 0, we assign them a group, otherwise we assign them a different group (i.e. we split into two groups). For each of these 2 groups we repeat the process, trying to split it further by recomputing Equation 6, looking at the signs of the leading eigenvector. If the leading eigenvalue is about 0 (or less), it is not a good division so that group of vertices is as maximally clustered as possible. We give these a unique label, and move to the next.

Answer 1

It's a bit late, but I have a very simple answer: there is nothing wrong with the code

I was simply missing the extra optimisation that Newman recommends in his paper. The results are perfectly in-line with what is reported in the paper without applying the extra optimisation step. I'll leave this up in case it helps someone implementing the algorithm.