Newman iteration foamalation

Hello,

I am trying to formulate an ILP formulation for the Newman algorithm step and would like to ask how can I improve it because it is very slow.

Newman algorithm's main idea is dividing a graph into two while maximizing modularity.

here are some definitions to understand Newman algorithm:

I got the help in this post to formulate this in Gurobi:

class Newman_ILP:
    def __init__(self, G, weight=None):
        """
        :param G: networkx graph
        :param nodes:
        """
        self.nodes_list = list(G.nodes)
        self.num_of_nodes = len(self.nodes_list)
        self.G = G
        self.weight = weight
        self.model = gp.Model("mip1")
        self.run()  # setting objective function and constraints and optimizing
        self.communities = self.get_communities()

    def run(self):
        self.set_objective()
        self.add_constraints()
        print("starting to optimize")
        self.model.optimize()

    """
    Objective Function: [sum_ij](q_ij * (y_ij + t_ijj))
    while: q_ij = a_ij - (d_i * d_j)/2m
    """

    @timeit
    def set_objective(self):
        m = self.G.size(weight=self.weight)
        adj = self.G.adj
        objective_function = 0
        for node_range_1 in range(self.num_of_nodes):
            i = self.nodes_list[node_range_1]
            globals()[f'x_{i}'] = self.model.addVar(vtype=GRB.BINARY, name=f'x_{i}')
            globals()[f't_{i}'] = self.model.addVar(vtype=GRB.BINARY, name=f't_{i}')
            for node_range_2 in range(node_range_1):  # i < j
                j = self.nodes_list[node_range_2]
                if dict(adj[i].items()).get(j) is not None:
                    a_ij = dict(adj[i].items())[j].get("weight", 1)
                else:
                    a_ij = 0
                q_ij = a_ij - (self.G.degree(i, weight=self.weight) * self.G.degree(j, weight=self.weight)) / (2 * m)

                globals()[f'y_{i}_{j}'] = self.model.addVar(vtype=GRB.BINARY, name=f'1_{i}_{j}')
                globals()[f't_{i}_{j}'] = self.model.addVar(vtype=GRB.BINARY, name=f'0_{i}_{j}')

                objective_function += (q_ij * (globals()[f'y_{i}_{j}'] + globals()[f't_{i}_{j}']))

        self.model.setObjective(objective_function, GRB.MAXIMIZE)

    @timeit
    def add_constraints(self):
        # Assumption: this function is called after set_objective() - which creates the variables
        for node_range_1 in range(self.num_of_nodes):
            for node_range_2 in range(node_range_1):  # j < i
                i = self.nodes_list[node_range_1]
                j = self.nodes_list[node_range_2]

                # y_ij ==1 iff x_j==x_i==1
                self.model.addGenConstrAnd(globals()[f'y_{i}_{j}'], [globals()[f'x_{i}'], globals()[f'x_{j}']],
                                           name="andconstr")
                # t_i ==1 - x_i
                self.model.addConstr(globals()[f't_{i}'] == 1 - globals()[f'x_{i}'])
                # t_j ==1 - x_j
                self.model.addConstr(globals()[f't_{j}'] == 1 - globals()[f'x_{j}'])
                # t_ij ==1 iff x_j==x_i==0
                self.model.addGenConstrAnd(globals()[f't_{i}_{j}'], [globals()[f't_{i}'], globals()[f't_{j}']],
                                           name="andconstr")

正式なコメント

Simranjit Kaur

Gurobi Staff

2025年11月16日 23:34

This post is more than three years old. Some information may not be up to date. For current information, please check the Gurobi Documentation or Knowledge Base. If you need more help, please create a new post in the community forum, or try Gurobot, our chatbot interface offering instant, expert-level support.

Maliheh Aramon

2022年06月18日 14:40

Hi Kim,

Your code is a bit complicated and I cannot read through it easily. As far as I understand, you are trying to implement the following mathematical model:

Notation
- \(G = (V, E)\): A graph with vertex set \(V\) and edge set \(E\)
- \(k_i\): degree of vertex \(i\)
- \(M\): \(2|E|\)
- \(A_{ij}\): 1 if edge \((i,j) \in E\)

Decision variable:
- \(s_i\): 1 if vertex \(i\) belongs to group 1 and -1 if it belongs to group 2

Objective function:
- \(\max Q = \frac{1}{2} \sum_{ij} (A_{ij} - \frac{k_ik_j}{M})(s_is_j+1)\)

A simple implementation is given below:

import gurobipy as gp
from gurobipy import GRB

import networkx as nx

if __name__ == "__main__":

     n = 20
     G = nx.erdos_renyi_graph(n, 0.5)
     A = nx.adjacency_matrix(G).toarray()
     k = [G.degree(i) foriinrange(n)]
     M = 2 * len(G.edges)

     model = gp.Model("modularity")
     x = model.addVars(n, vtype=GRB.BINARY, name="x")
     s = model.addVars(n, lb=-1, ub=1, vtype=GRB.INTEGER, name="s")

     # Define s_i=2x_i -1 such that if x_i=1, then s_i=1 and
     # if x_i=0, then s_i=-1
     model.addConstrs(
          (s[i] == 2 * x[i] - 1 for i in range(n)), name="connect_x_s_constr"
     )
 
     # Define the objective function
     Q = gp.QuadExpr()
     for i in range(n):
         for j in range(n):
             Q += (A[i, j] - k[i] * k[j] / M) * (s[i] * s[j] + 1)  
      Q *= 0.5
      model.setObjective(Q, sense=GRB.MAXIMIZE)
      
      model.optimize()

Best regards,

Maliheh

Newman iteration foamalation

コメント

お探しのものを見つけられませんでしたか？