Skip to main content

Newman iteration foamalation

Answered
Kim Kestenboum
  • Gurobi-versary
  • First Comment
  • First Question

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")
0

Comments

2 comments

Date Votes
  • Official comment
    Simranjit Kaur
    • Gurobi Staff
    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
    • Gurobi Staff

    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

    1

Post is closed for comments.

Didn't find what you were looking for?

New post