New Advances on Analysis and Optimization Problems in Networks
Abstract
En esta tesis doctoral, se desarrollan nuevos avances en problemas de clasificación de los elementos de una red, con el fin de analizar en profundidad y comprender mejor este tipo de estructuras y sus variantes más complejas. Estos problemas se abordan siempre desde el punto de vista de la optimización, a través de funciones representativas o medidas de calidad de la clasificación de los nodos de la red realizadas. De esta forma, podemos apoyarnos en la gran variedad de herramientas que nos proporcionan campos como la Investigación Operativa y la Programación Matemática para resolver este tipo de problemas. El documento se divide en dos partes. La primera parte consta de cinco capítulos que introducen, argumentan y presentan ciertos aspectos claves para comprender el contenido teórico y científico de la tesis que se desarrolla en profundidad en la segunda parte de esta. El Capítulo 1 introduce el marco general y los conceptos básicos necesarios para comprender el contexto en el que se ubica esta tesis, así como motivarla. En el Capítulo 2, destacamos los objetivos y puntos fuertes que se han perseguido durante la realización de la tesis. El Capítulo 3 presenta los resultados obtenidos y objetivos que han sido alcanzados. La importancia y consecuencias de dichos resultados se discuten en el Capítulo 4. Finalmente, el Capítulo 5 presenta las conclusiones generales del trabajo y posibles futuras líneas de investigación que pueden surgir a raíz de los resultados obtenidos en esta tesis. La parte dos está compuesta por cuatro capítulos, Capítulos 6-9, que consisten en cuatro trabajos de investigación independientes y autocontenidos realizados durante el desarrollo de la tesis. Además, tres de los cuatros capítulos equivalen a artículos de investigación publicados durante el proceso de tesis en medios relacionados con el campo científico de esta.
Full text
PHD DISSERTATION NEW ADVANCES ON ANALYSIS AND OPTIMIZATION PROBLEMS IN NETWORKS. Francisco Temprano García Supervisor: Prof. Dr. Justo Puerto Universidad de Sevilla Instituto de Matemáticas de la Universidad de Sevilla (IMUS) Programa de Doctorado en Matemáticas Sevilla XX.XX.2025
A mis padres y mi hermana. III
Agradecimientos En primer lugar, me gustaría agradecer a mi director de tesis Justo Puerto Albandoz por ser el principal causante de la realización de esta tesis. Gracias por brindarme esta oportunidad y mostrarme el camino de la carrera investigadora. Si no es por Justo, yo no habría hecho esta tesis ni me dedicaría a esto, él es el verdadero y único artíce. Gracias también por entender mi forma de trabajar y dedicación, yo sé mas que nadie lo complicado que es colaborar y trabajar con una persona como yo. También, me gustaría agradecer a todos mis coautores que han aportado su sabiduría y conocimientos a este trabajo, y de los que he aprendido mucho al trabajar con ellos. Gracias a Antonio Manuel Rodríguez Chía de la Universidad de Cádiz, Stefano Benati de la Universidad de Trento y Diego Ponce de la Universidad de Sevilla por vuestra ayuda y dedicación, así como la compañía y momentos juntos. Además, me gustaría agradecer tanto a Stefano como a Lina Mallozzi de la Universidad de Nápoles Federico II por permitirme realizar mis estancias investigadores del doctorado en Trento y Nápoles, y hacerme sentir como en casa a pesar de estar en otro país, aunque maravilloso, como es Italia. Dejando a un lado a las personas que han contribuido de forma directa en esta tesis, hay muchas más personas que me han acompañado durante estos años y han contribuido de otra forma. Especial dedicación a mis dos apoyos más importante en el mundo de las matemáticas, Nacho y Alberto. Hace prácticamente 10 años que comenzamos juntos el Doble Grado de Matemáticas y Estadística. Durante esos 10 años no nos hemos separado y hemos tomado casi las mismas decisiones y caminos. Para mí es un sueño que los tres hayamos decidido hacer el doctorado, y no sé si hubiera sido capaz de hacerlo solo. También me gustaría agradecer a aquellas personas que desde una posición más cercana me han servido de guías y modelos a seguir. Tener a personas como Carlos y Moisés cerca me hizo aprender rápido lo que es el doctorado y la investigación. Ellos me enseñaron tanto lo bueno como lo malo de la vida académica. Tener a dos de los mejores investigadores de mi área y personas más inteligentes que he llegado a conocer en tu mismo despacho es una gran ventaja. Incluso a día de hoy, he tenido de modelos sus dos tesis doctorales para la redacción de la mía. Moi, gracias por tanto, te echo de menos. No me puedo olvidar de todos mis compañeros y personas que han pasado por el Instituto de Matemáticas de la Universidad de Sevilla (IMUS). Los días, comidas, cafes, risas y momentos compartidos con todas las personas que he conocido aquí son incontables. Son muchísimas las personas que considero mis amigos y he compartido momentos dentro y fuera del IMUS, pero querría agradecer el apoyo y compañía a Juan, Carlos Madrid, V
VI Carlos Pelirrojo, Nati, Cristina, Reme, Tom, Bandera, Nuria, Jasone, Gabri, Manu... Sé que me dejo a mucha gente, pero muchas gracias de verdad a todos. Durante el doctorado he podido disfrutar de muchos viajes y congresos donde he conocido a mucha gente, de las que me llevo a muchos amigos. Esto me ha servido para conocer a grandes investigadores con los que comparto una misma pasión y amor por las áreas de Investigación Operativa, Programación Matemática y Optimización. Gracias a ellos siento que pertenezco a una gran familia y cada vez que voy a algún congreso es una reunión con personas que respeto y aprecio. A pesar de no tener ningún tipo de relación con el mundo académico, me quiero acordar y agradecer de todo aquel que me ha acompañado durante mi doctorado. Entre ellos se encuentran todos mis amigos que llevan conmigo prácticamente toda la vida, y que a día de hoy siguen siendo pilares fundamentales. Rau, Luis, Marce, Morci y todos mis amigos de Mairena o cercanía, todo es más fácil, más divertido y mejor gracias a vosotros. Gracias también a los amigos que hice durante mi época de estudiante universitario y con los que sigo manteniendo la amistad. Por último pero no menos importante, quiero agradecer a toda mi familia a la que amo y de la que me siento tan orgulloso. Espero que nunca me falten, y seguir sintiendo su amor incondicional toda la vida. También acordarme de aquellos familiares que ya no están, el orgullo que siento de ser nieto de mis abuelas es eterno. Mama, Papa, Isa, familia, amigos y todas las personas que quiero, todo lo que hago y soy es por y gracias a vosotros.
Resumen En esta tesis doctoral, se desarrollan nuevos avances en problemas de clasicación de los elementos de una red, con el n de analizar en profundidad y comprender mejor este tipo de estructuras y sus variantes más complejas. Estos problemas se abordan siempre desde el punto de vista de la optimización, a través de funciones representativas o medidas de calidad de la clasicación de los nodos de la red realizadas. De esta forma, podemos apoyarnos en la gran variedad de herramientas que nos proporcionan campos como la Investigación Operativa y la Programación Matemática para resolver este tipo de problemas. El documento se divide en dos partes. La primera parte consta de cinco capítulos que introducen, argumentan y presentan ciertos aspectos claves para comprender el contenido teórico y cientíco de la tesis que se desarrolla en profundidad en la segunda parte de esta. El Capítulo 1 introduce el marco general y los conceptos básicos necesarios para comprender el contexto en el que se ubica esta tesis, así como motivarla. En el Capítulo 2, destacamos los objetivos y puntos fuertes que se han perseguido durante la realización de la tesis. El Capítulo 3 presenta los resultados obtenidos y objetivos que han sido alcanzados. La importancia y consecuencias de dichos resultados se discuten en el Capítulo 4. Finalmente, el Capítulo 5 presenta las conclusiones generales del trabajo y posibles futuras líneas de investigación que pueden surgir a raíz de los resultados obtenidos en esta tesis. La parte dos está compuesta por cuatro capítulos, Capítulos 6-9, que consisten en cuatro trabajos de investigación independientes y autocontenidos realizados durante el desarrollo de la tesis. Además, tres de los cuatros capítulos equivalen a artículos de investigación publicados durante el proceso de tesis en medios relacionados con el campo cientíco de esta. En los Capítulos 6 y 7, abordamos desde dos puntos de vista totalmente distintos la extensión del problema clásico de detección de comunidades donde el solapamiento entre los distintos subgrupos en los que clasicamos los nodos de una red es posible. El objetivo principal de ambos trabajos es proporcionar una metodología, basada en Programación Matemática, que evalúe la bondad de las medidas de calidad propuestas en la literatura y que, además, nos permita proponer variantes que corrijan las desventajas de estas. El Capítulo 6 se enfoca en aplicar dicha metodología a un trabajo, de gran peso en la literatura de estos problemas, que introduce el concepto de anidad de un elemento a su comunidad, que se puede entender como la probabilidad de pertenencia. Gracias a nuestra metodología, llegamos a la conclusión de que la medida de calidad propuesta no es realmente una función representativa de las clasicaciones, y basándonos en el concepto de anidad y modularidad proponemos una medida de calidad alternativa. A través de VII
VIII experimentos computacionales, validamos las ventajas y buenas propiedades de nuestra nueva función objetivo. Además, diseñamos varios algoritmos heurísticos para este problema de optimización, cuya eciencia y comportamiento es evaluado mediante distintos experimentos sobre algunas redes conocidas en la literatura de este campo. Por otro lado, el Capítulo 7 aplica nuestra metodología sobre otro trabajo bastante conocido que se olvida del concepto de anidad pero se apoya en muchos aspectos, propiedades y resultados de la Teoría de Juegos Cooperativos. En este caso, proponen una medida de calidad que consiste en maximizar la suma de los valores de Shapley de las coaliciones detectadas del juego cooperativo, que, además, deberán garantizar una conocida propiedad llamada estabilidad. Gracias a los modelos de optimización, podemos evaluar la bondad de dicha medida de calidad, y proponer nuevas funciones que mejoren sus resultados. En este trabajo, basándonos de nuevo en la modularidad, proponemos una nueva medida de calidad que es valorada mediante extensos experimentos computacionales. Además, un heurístico es diseñado para poder obtener soluciones de buena calidad en instancias de tamaños superiores. El rendimiento y eciencia de todos los métodos que se proponen son comparados sobre instancias de redes muy conocidas en la literatura. Dejando a un lado los problemas de detección de comunidades solapadas, el Capítulo 8 se concentra en un problema de optimización bastante conocido en la literatura y en el área de problemas de agrupamiento: el problema de los k -cortes normalizados mínimo. El objetivo de dicho problema es encontrar la k -partición que minimice la densidad de aristas entre los nodos que son clasicados en subgrupos distintos. A pesar de que este problema ha sido estudiado por muchos investigadores durante más de vente años, nadie ha propuesto una formulación o modelo capaz de resolver de forma exacta este problema de optimización. En este trabajo, presentamos una extensa variedad de métodos capaces de encontrar el óptimo exacto. Estos métodos consistirán en formulaciones compactas, que pertenecen a la familia de problemas de programación lineal entera-mixta, y un algoritmo basado en generación de columnas que también es capaz de resolver de forma exacta el problema. Todas las familias de métodos son comparadas de forma justa mediante extensos experimentos computacionales sobre instancias generadas aleatoriamente. Finalmente, proponemos un algoritmo heurístico capaz de superar las desventajas de escalabilidad de los métodos exactos y que puede ser aplicado como algoritmo de segmentación de imágenes. Este algoritmo es aplicado sobre una serie de imágenes que aparecen en la literatura del problema y de las que se quiere proporcionar una división de forma que cada parte tiene un signicado propio. Por último, el Capítulo 9 estudia los problemas de detección de comunidades para el caso especial de hipergrafos. Un hipergrafo es una estructura más compleja que una red simple donde los patrones de conexión pueden consistir en cualquier subconjunto de nodos de cualquier tamaño, a diferencia que en un grafo simple donde las conexiones son siempre entre pares de nodos. Aunque los hipergrafos son estructuras que han sido analizadas y son utilizadas para representar una extensa variedad de tipos de datos, la gran mayoría de modelos y métodos de agrupamiento sobre ellos se basan en transformarlos en grafos simples u otras estructuras sobre las que aplicar otros métodos ya conocidos.
IX Es por esto, que proponemos un nuevo modelo que trata el hipergrafo desde su estructura original mediante la extensión de ciertos conceptos ya existentes para el caso simple. En concreto, extendemos el concepto de hiperarista interna, su grado de internalidad y modelo de hipergrafos aleatorios. Los ingredientes anteriores nos permiten presentar una nueva denición de la función de modularidad en hipergrafos. Este nuevo problema de optimización es modelado mediante una formulación lineal entera mixta que nos permite resolverlo de forma exacta. La bondad y el rendimiento del método son validados mediante experimentos computacionales sobre hipergrafos generados aleatoriamente. Con el objetivo de mostrar la gran aplicabilidad del problema, hemos clasicado las posibles respuestas de dos datos de encuestas pertenecientes al Eurobarómetro y al Barómetro del Centro de Investigaciones Sociológicas (CIS) de España.
Chapter 1 Introduction 1
3 A network is a complex structure constituted by two key components: a set of elements called nodes, that can be understood as a population; a set of connections or edges between nodes that dene the dierent existing relationships and patterns between the items of the node set. While the classic node set denition is universally accepted and admitted, the nature of the connection set can dier in many ways depending on the type of structure that one wants to dene and the benets that one wants to obtain. These structures have been used to represent huge datasets for a wide variety of subareas such as biology, computer science, location, transportation, sociology, economics, and many more (see Newman (2010)). Once data is represented as a network, it is possible to take advantage of it and obtain conclusions about it. Here, Operations Research (OR) and Mathematical Programming (MP) play a key role as decision tools for the dierent problems that arise. These techniques have been widely applied in well-known discrete optimization problems characterized by networks, such as facility location (Nickel and Puerto, 2009), hub location (Saldanha-da Gama and Wang, 2024), logistic (Fischetti et al., 2006), supply chain management (Stevens, 1989) or game theory (Owen, 1982). Nowadays, there are countless real-world challenges on which applying techniques, models and algorithms, as those mentioned above, can help us to provide good quality solutions to them. Among these problems, we can highlight cost savings, improved quality of services and better organization of resources and population structures, always through large amounts of data that serve as a basis for facing these challenges. In conclusion, the use of data to make optimal or good quality decisions in our favor is becoming more and more important. For all these reasons, network analysis and optimization play a fundamental role in making these kind of decisions. These subjects provide the mathematical framework and methodology that allow the resolution of the dierent considered problems. In simple terms, an optimization problem consists of the maximization or minimization of an objective function that depends on a set of decisions variables that must belong to a restricted subset. This restricted subset is called feasible decision subset, since it includes all the possible decisions of which we want to nd the optimal one. This is where Mathematical Programming and Operational Research come in to be able to model and formulate an optimization problem mathematically through its objective function and its set of constraints that dene the denitive feasible subset. Hence, all the known tools and machinery can be used to solve any optimization problem formulated by Mathematical Programming. We refer above to the state-of-the-art algorithms and solvers implemented in any of the programming languages we can think as C, C++, C#, Java, Visual Basics and Python, that are able to provide decent resolutions for instances up to considerable sizes. This last assessment is due to the complexity limitations of the problems that are well known in our reserch area and have not been possible to overcome until now. However, this does not mean that there are no ways to improve the performance of the resolution methods and formulations already proposed in literature. The in-depth study and analysis of a problem formulation or structure can reveal specic characteristics and properties that can help to
4 Chapter 1. Introduction improve the results of the already known methods and reach new heights of complexity that were not possible before. This thesis applies Mathematical Programming models and algorithms to networks and other data represented over a graph with the aim of providing a classication of the elements that belong to these structures. These individuals or members can be organized in many possible ways, we will take into account from the most classic denition to the most complex structures and extensions. In general, all the studied problems consist of providing subgroups of the element set, called communities or clusters, such that members from the same community are considered similar or with similar features, and members that belong to dierent communities have dierent characteristics. This kind of problems are commonly called Community Detection problems which receive wider attention thanks to the work of Girvan and Newman (2002), although clustering problems go back longer with the sociology, zoology, economics, and political sciences applications of Grötschel and Wakabayashi (1989); Grötschel and Wakabayashi (1990); Wasserman and Faust (1994). A key point of the Community Detection problems is the way similarity between nodes is measure. In the most classic denitions, the similarity between nodes from the same community is measure by the internal edge density of these subgroups. However, it is possible to nd alternatives on which the aim is to classify in the same community nodes with similar connection patterns. Another assumption that will be repeated throughout the document will be that each element must belong to at least one of the communities, i.e., the subgroup set is a subdivision of the network elements. The main purpose of the thesis is to design formulations, methods and algorithms capable of solving the above combinatorial problem through community structures that organize properly the elements of a network based on the connections and relationships of the network and other factors that provide us with information of our data. This will be possible thanks to in-depth study of the properties and structures of the problems that we have addressed through a Mathematical Programming and Operations Research point of view. In this rst section, the necessary background to comprehend properly and motivate this family of classic combinatorial problems is presented. This includes from the most general theory of Mathematical Programming to those classical problems that are more related to the Community Detection eld. In addition, we will introduce a basic knowledge of complexity theory due to the great relationship that it has with this kind of optimization problems. Section 1.1 focuses on the most classic Mathematical Programming theory and explain its most basic concepts. Next, Section 1.2 presents some of the most known Community Detection methods which also have a great inuence and impact on this work. Finally, Section 1.3 gives an idea of Computational Complexity Theory and its relationship with optimization.
1.1. Optimization problem 5 1.1 Optimization problem An optimization problem is dened as the problem of nding the best solution from all feasible solutions with respect to an objective function. The standard form of an optimization problem is: minimize xf(x) (OP) subject to gi(x)≤0,∀i= 1, . . . , m, (1.1a) hj(x) = 0,∀j= 1, . . . , m, (1.1b) x∈X⊆Rn, (1.1c) where f:Rn→R is the objective function to be minimized over the n -variable vector x . The feasible set is dened by the constraint functions gi and hj , and the domain X⊆Rn that can be both continuous and discrete. This is the most general and primitive form of an optimization problem, since any problem can be expressed as (OP). Although this denition includes an immense variety of problems, there are properties and methods to solve them, such as the well-known Karush-Kuhn-Tucker (KKT) conditions (see Karush (1939); Kuhn and Tucker (1951)). Within all the classes of optimization problems, we can higlight Linear Programming (Dantzig et al., 1954), Quadratic Programming (Markowitz, 1952), Second Order Cone Programming, Semidenite Programming and Conic Programming (Nesterov and Nemirovski, 1992; Nesterov and Nemirovskii, 1994; Alizadeh, 1995). Most of the models developed in the work belong to the class of Linear Programming, and more specically to the Mixed-Integer Linear Programming (MILP). This is due to the combinatorial nature of Community Detection problems. In the following, this family of problems and its most important properties are introduced. 1.1.1 Linear Programming A Linear Programming (LP) problem is an optimization problem that can be expressed in standard form as: minimize xctx (LP) subject to Ax =b,, (1.2a) x≥0, (1.2b) where c∈Rn and b∈Rm are given vectors, A∈Rm×n is a given matrix and (A, b)∈ Rm×(n+1) is a full rank matrix. The problem (LP) can be also dened as nding the solution x in the polyhedron Ax =b that minimizes the linear expression ctx . Firstly, we can highlight that the solution of (LP) must be an extreme point of the feasible set Ax =b , thanks to the linearity of the problem. Dantzig et al. (1954) presented the well-known Simplex algorithm to solve exactly any LP problem, based on going through the extreme points of the polyhedron, improving the objective function until reaching the
6 Chapter 1. Introduction minimum. Moreover, it is known that every extreme point in the polyhedron is univocally associated with a subset of variables of size m , so it only depends on this subset and can be obtained by solving the equation system Ax =b with this reduced subset of variables. This means that a Linear Programming problem can really be solved only with the subset of variables associated with the optimal extreme point. The variable subset associated to an extreme point is usually called basic variables. However, there is not always a solution to these problems, since there may be two other cases: a LP is infeasible if the feasible set {x∈Rn:Ax =b} is empty; a LP is unbounded if the objective function could be reduced innitely within the feasible set. Moreover, for any LP problem there is an equivalent LP problem called Dual problem. In the rest of the section, we will refer to (LP) as Primal problem. The dual problem of (LP) is dened as follows: maximize ubtu (D) subject to Atu≤c, , (1.3a) u∈Rm, (1.3b) The equivalence and properties of both Primal and Dual problems are key to the solving algorithms of this family of problems. The importance of duality in Linear Programming is summarised in the following theorem (Wolfe, 1961): Theorem 1 Given a primal Linear Programming problem (LP) , and its dual (D) , they verify: 1. Weak Duality : For any feasible solution x of (LP) , and any feasible solution u of (D) , then: ctx≥btu. 2. Strong Duality : The optimal values of (LP) and (D) are the same, i.e., if ¯x and ¯u are the optimal solution of (LP) and (D) , respectively, then: ct¯x=bt¯u. 3. Exactly one of the following three cases occurs: (a) Both (LP) and (D) are infeasible. (b) One of (LP) and (D) is unbounded, and the other is infeasible. (c) Both (LP) and (D) have optimal solutions, and their optimal values are identical. The above theorem and equivalences between (LP) and (D) are only valid when the variables of (LP) are dened over a continuous domain. In the following section, we will show the classical clustering formulations and problems that belong to the MixedInteger Linear Programming family of problems. These problems belong to the most
1.2. Community detection problems 7 complex type of problems to solve, the NP-hard class, (see Cook (1971); Karp (1972)), while the continuous (LP) problems are known to be solvable in polynomial time thanks to the algortihm proposed by Khachiyan (1980). However, even with MILP problems, it is possible to take advantages of the above theorem and the duality properties, so that the performance of the existing algorithms is improved. Methods based on the above duality theorem, as the Lagrangian Relaxation (Everett, 1963) and Benders Decomposition (Benders, 1962), works well and are really useful in the resolution of MILP problems. Returning to the concept of basic variables, there are a large number of exact algorithms that try to reduce the dimension of the variable set in order to make the problem easier but preserving the optimality of the solution obtained. In the best case, the algorithm solves the problem only with the subset of basic variables of the optimal solution, and, in the worst case, the algorithm needs to solve the problem with the complete variable set. Due to the fact that the optimal basic variables subset is unknown a priori, every time the problem is solved for a reduced variable subset it is necessary to check whether the obtained solution is optimal or not. Formally, the following proposition can be declared: Proposition 1 Given an extreme point x∗ of the polyhedron {x≥0 : Ax =b} represented by the basic variables xB and the non-basic variables xN , then it is known that x∗ B= B−1b≥0 and x∗ N= 0 , where B∈Rm×m is a nonsingular submatrix consisting of the columns of A associated with basic variables xB . In addition, we can arm that the extreme point x∗ is the optimal solution of (LP) if and only if: c−cBB−1A≥0, (1.4) where cB are the values of c associated to the basic variables xB . The left hand side vector of (1.4) is called reduced costs and if it contains at least one negative value, then there is a dierent extreme point with a strictly lower objective value. Thanks to the above proposition, it is possible to verify the optimality of the solution obtained by a reduced subset of variables. Algorithms as Column Generation for LP problems or Branch-and-Price for MILP problems are based on this key proposition, so they know when the optimal solution is reached, or if, on the other hand, it is necessary to add new variables to the pool of the actual considered variables. 1.2 Community detection problems The classic Community Detection problems consist of provinding a partition of the node set of a graph or network, in such a way that nodes from the same subgroup have similar features and nodes from dierent clusters do not. Thus, this partition, called Community Structure, organizes properly this population that is represented as a network, and allows us to better understand the data and obtain correct conclusions. Although all Community Detection problems have the same aim and structure, there are factors and details that make them dierent from each other. A key factor is the structure of the algorithm, for example, Hierarchical Clustering, as the proposed by Sibson
8 Chapter 1. Introduction (1973), is a well-known family of algorithms whose objective is to design a dendrogram. A dendrogram is a diagram representing a tree that starts from the trivial partition, formed by the total set of nodes, and makes bipartitions that generate new partitions until reaching the partition formed by single node subgroups. So, instead of obtaining a single partition, the dendrogram provides a family of possible community structures. Within all classes of clustering algorithms we will focus on those that are based on designing an objective function to optimize. It is really common to develop a quality measure as a representative value of the community structure, so that the optimal solution is equivalent to the partition that best classies the node set into communities. Once we turn the Community Detection problem into a Optimization Problem, we can apply all the tools and algorithms provided by the Mathematical Programming to solve it exactly or obtain a good approximation. In the following, we present three Community Detection problems that belong to this class of Clustering algorithm and are among the best-known and most used in literature. In addition, most models, algorithms and theory developed in this document are based on these three problems. 1.2.1 Clique Partition Problem Let V={1, . . . , n} be a node set, and D= (dij ∈R)i,j∈V i<j be a weight set that measures the disimilarity or distance between nodes for any possible pair of the set V . The Clique Partition Problem (CPP) consists of nding the partition P={C1, . . . , Cq} of V that minimizes the objective function: q X s=1 X i,j∈Cs i<j dij. (1.5) Due to the problem denition, we can interpret that two nodes i and j are similar if dij ≤0 , and they are not if dij >0 . So, the disimilarity weight set D can take both positive and negative values. Once the CPP is presented, Grötschel and Wakabayashi (1989) proposed the following Mixed-Integer Linear Programming formulation to solve exactly this problem: min n−1 X i=1 n X j=i+1 dijxij (1.6a) s.t. xij +xjk −xik ≤1, i, j, k = 1, . . . , n, i < j < k, (1.6b) xij −xjk +xik ≤1, i, j, k = 1, . . . , n, i < j < k, (1.6c) −xij +xjk +xik ≤1, i, j, k = 1, . . . , n, i < j < k, (1.6d) xij ∈ {0,1}, i, j = 1, . . . , n, i < j, (1.6e) where xij is the binary variable that takes value 1 if nodes i and j are in the same community, and takes value 0 otherwise. Expression (1.5) is minimized through the objective function (1.6a). In addition, binary x -variables must guarantee constraints (1.6b)-(1.6d) to be well-dened, since they must satisfy the transitive property. The transitive property
1.2. Community detection problems 9 means that if nodes i and j are in the same community, and nodes j and k are in the same community, then nodes i and k have to be also classify in the same subgroup. In Grötschel and Wakabayashi (1989); Grötschel and Wakabayashi (1990), it is also deeply analyzed the CPP polytope described in (1.6). They found some facets and polyhedric properties in order to described the convex hull of the feasible solutions of this problem in more detail. In Grötschel and Wakabayashi (1989), we can obserse a large number of real databases in which the CPP is applied, showing the great applicability of this clustering algorithm. It is possible to nd classications of animals and plants species (Biology), classications of people populations (Sociology), classications of vehicles and other type of structures (Engineering), classications of computers and other type of machines (Computer Science), and classications of Nations or other political systems (Political Science). 1.2.2 Modularity One of the most well-known quality measure in network and clustering analysis that measures the goodness of our community structure is the so-called modularity. This expression was introduced by Newman and Girvan (2004) in order to provide a meaningful classi- cation for the items/elements of a graph or network when the relationships between these elements are dened by the edges of the graph. Although the modularity was originally proposed for networks, it can be applied in a wide variety of data structures with countless advantages. This is due to the fact that there are many techniques to turn all these data structures into networks, so we can apply this modularity concept. Given an undirected graph G= (V, E) , where V={1, . . . , n} is the node set, and E is a set of node pairs called edges that represent the connections and relationships between nodes. The aim of modularity is to provide a function that measures the internal edge density of the communities we detect. Newman and Girvan (2004) considered that nding a good community structure is equivalent to arm that the edge set are not randomly distributed through the whole node set, since the probability of connecting nodes from the same community is larger. So, modularity compares the actual internal number of edges between nodes from the same community of our detected structure with its expected value in random graph that follows the random distribution proposed in Newman (2010). These random graphs are built by splitting each edge into two stubs, each one adjacent to one of the nodes of the pair, and then randomly joining these stubs pair by pair. So, this random graph guarantees that it has the same number of edges and their nodes have the same adjacent degree than the original graph G . In addition, we can assume that each random graph is uniquely dened by the matching that joins pair by pair the stub set. Once the random distribution is dened, it is possible to exactly calculated the expected number of edges between any pair of nodes. Formally, for any nodes i and j let Aij be the number of edges in E between i and j . Thus, the adjacent degree of a node i , dened as the number of edges connected with i , can be expressed as ki=P j∈V Aij . Now, given
17 This chapter declares the main objectives of this thesis. As it has been commented in the introduction, the aim of this work is the use of Mathematical Programming to extend the applicability of Community Detection methods to new interesting and realistic scenarios, as well as develop theorical results on the models already proposed in the literature. These objectives are summarized in the following bullet points: To introduce new extensions or variations of the classic Community Detection problems with the aim of achieving more complex situations and structures in which there is not universal agreement or no proposals have been made. This will always be done in search of realistic scenarios with great real applicability. To provide a methodology to solve clustering problems from a Mathematical Programming point of view. We will rely on the tools and algorithms, that are capable of exactly solving optimization models, to address these problems. This is due to a signicant lack of exact resolution methods in the literature of this nature. To develop theorical results and advances in the Operations Research knowledge of these problems through an in-depth study of their properties and structures that can help us to simplify them and speed up their resolution. In addition, we will try to design heuristic algorithms to complement our exact approach and provide good quality approximations to the optimal solution in a signicantly shorter time. To test and compare computantionally the proposed methods with the already existing ones by means of their performances, eciencies and structures. Additionally, we will apply our methods to real data and instances, in order to show the great applicability and impact of clustering and Community Detection problems in the society and real scenarios.
Chapter 3 Results 19
21 The results obtained during the development of this thesis are exposed in this chapter. In addition, these results have resulted in three publications: Benati et al. (2022), Benati et al. (2023) and Ponce et al. (2024); and a preprint: Benati et al. (2024); which correspond to Chapters 6, 7, 8 and 9, respectively. Obviously, they are strongly related with the objectives presented in the above section, and can be summarised as follows: We have introduced new adaptations of some problems in the Community Detection eld with a great interest and applications. Specically, we have tried to solve the lack of universal agreement on Overlapping Clustering through Mathematical Programming. This is an extension of the classic problems where overlap between communities is allowed, which, in certain situations, are more realistic than the usual disjoint structures. The methods presented in literature do not seem to be accurate enough, that is why, we have tested their goodness and provided the evidences that our new proposals can x their biases. In addition, we have proposed new methodologies to solve clustering problems over complex network structures that were not considered before, such as hypergraphs in which the size of their edges and connections could be larger than two. Mixed Integer Linear Programming models have been developed to exactly solve the above extensions and some classic problems in literature for which there were no exact resolution methods. The latter refers to the MkWNCP, for which we have presented some MILP formulations and a Column Generation algorithm capable of solving it. Moreover, for every problem considered in this thesis, some heuristic algorithms have also been designed due to the complexity and scalability disadvantages of the exact approaches. Heuristics have allowed us to carry out realistic applications of great impact. Hand in hand with the above achievements and models, we have obtained some theoretical results with great inuence in the elds of Mathematical Programming and Operations Research. This may help in future works about Community Detection or other optimization problems. Among them we can highlight NP-hardness, polyhedric results and other interesting properties about the structure of these problems. All the above presented material is contrasted through an extensive set of computational experiments that have allowed us to compare and test the performance of the large number of algorithms and methods proposed in this document. We have also compared ourselves with the existing literature to show the contribution of this thesis. These comparisons must be made fairly, so, the generation and selection of the instances on which the experiments have been applied are also key points.
Chapter 4 Discussion of the results 23
25 In this chapter, the contribution and novelty of the obtained results, that are presented in the above section and explained in detail in Chapters 6-9, are going to be argued and supported. Chapter 6 is the rst of the works that will serve to propose and show the necessity for a new family of methodologies in the state-of-the-art Community Detection eld based on Optimization and Mathematical Programming. Given a variation or related Community Detection problem based on optimizing a designed quality and representative measure, this methodology consists of using the tools and models provided by Mathematical Programming to achieve two clear objectives: to test the goodness, eectiveness and performances of the algorithms and quality measures already proposed in literature, since we are able to obtain their optimal solutions; to develop models, formulations and methods capable of exactly solving these problems as ecient as possible. Specically, in this paper we focus in one of the most cited and known quality measures for overlapping community structures in literature. We talk about the fuzzy modularity function that was originally proposed by Zhang et al. (2007) which introduced a new concept of membership of an element to a community and extended the modularity measure to overlapping structures through this new ingredient. However, no attempt was made to provide an exact algorithm for the optimization problem, and they just presented a heuristic algorithm to supposedly approximate the optimal solution. Through our work, we formalize the optimization problem as a Mixed Integer Quadratic Programming model and show the biases of the objective function proposed by Zhang et al. (2007). After that, we provide a dierent fuzzy modularity function based on similar resources that x the misbehaviour of the original measure and can be modeled as a MILP. We prove the superior eciency of our function by comparing both exact models over some well-known data sets of great impact in the literature of these problems. Finally, three heuristics algorithms based on local search are proposed and applied over the same data sets. Those experiments show their good performances and that properly designed heuristics are able to solve the scalability disadvantages of the exact model through really good quality approximations. Chapter 7 reapplies the methodology proposed in the above paragraph to another wellknown overlapping clustering method. In this case, we refer to the algorithm proposed by Jonnalagadda and Kuppusamy (2018). In this work, a cooperative game, called weighted graph community game, is dened so that each stable coalition of the game can be seen as a possible community of the overlapping structure. They also propose to maximize the sum of the communities Shapley value, one of the most representative index in cooperative game theory (Shapley, 1951). Again, thanks to formalize the optimization problem as MILP model, we point out the misbehaviour and biases of the algorithm proposed by Jonnalagadda and Kuppusamy (2018). So, we try to apply the modularity denition to the similarity weights that appear on the weighted graph community game. These weights are dened in a complex way, so the exact modularity expresion can not be calculated trivially. Some probability results and expressions are exactly calculated, such as the adjacency probability, using the properties and combinatorial structure of the random graph generation that was presented in Chapter 1. Once this new quality version is
32 Chapter 5. Conclusions the random instances generator presented in Lancichinetti et al. (2008) that will allow anyone to test and compare fairly any kind of algorithm over these instances. Thanks to extensive computational experiments, we assess the strong and weak points of each one of the methods proposed in this work. As future work, we can point out that stability is a really interesting concept in clustering and cooperative game elds that is poorly analyzed in literature and increases signicantly the problem complexity. In Chapter 8, we propose multiple solution methods to solve the MkWNCP which has been analyzed in the last 20 years and whose advantages and benets are well-known. However, no attempts have been made to solve exactly this problem before. We propose for the rst time a MILP reformulation of the original non-linear normalized cut function. Furthermore, a Branch-and-Price, based on Column Generation and Set Partition formulation, is developed to show the eciency of algorithms of this nature. The performances of all these methods are tested with randomly generated instances which allows us to guess the best models. In addition, a matheuristic, that use partially the above exact methods, is apply as image segmentation algorithm to the same images of the paper of Shi and Malik (2000). Thanks to these experiments, we can conclude that our matheuristic behaves properly and outperforms other well-known heuristics, as the one proposed in Shi and Malik (2000). As a possibe improvement, we can comment that Column Generation procedure can be accelerated by some tricks and properties of the problem structure. Chapter 9 proposes a new denition of hypergraph modularity as a function capable of measuring the goodness of communities dened on hypergraphs. Most of the methods in literature consist of projecting hypergraphs to simple graphs to which apply all the already known classic algorithms for simple networks. Nonetheless, we think this is inconsistent because of several reasons, and our contribution deals directly with the original hypergraph. We develop a new random hypergraph model and a new concept of internality that allow us to develop a hypergraph modularity function and formulate it as a MILP model. We compare both approaches, the one that projects the hypergraph and our proposed modularity function. We develop also a novel hypergraph benchmark, that allow us to conclude that our method outperforms the projection approach over these instances randomly generated. Finally, we show a peculiar and interesting application over some survey data, in which the answer list can be modeled as hypergraphs. The results point out that modularity can cluster them more accurately than other standard statistical techniques. As a possible improvement, we have not focused on the designing of any heuristic to optimize the hypergraph modularity. So, the need to implement heuristic algorithms for large size hypergraph is clear, due to its resolution complexity.
Bibliography Agarwal, G. and Kempe, D. (2007). Modularity-maximizing graph communities via mathematical programming. Physics of Condensed Matter , 66. Alizadeh, F. (1995). Interior Point Methods in Semidenite Programming with Applications to Combinatorial Optimization. SIAM Journal on Optimization , 5(1):1351. Arenas, A., Duch, J., Fernández, A., and Gómez, S. (2007). Size reduction of complex networks preserving modularity. New Journal of Physics , 9(6):176. Benati, S., Puerto, J., Rodríguez-Chía, A. M., and Temprano, F. (2022). A mathematical programming approach to overlapping community detection. Physica A: Statistical Mechanics and its Applications , 602:127628. Benati, S., Puerto, J., Rodríguez-Chía, A. M., and Temprano, F. (2023). Overlapping communities detection through weighted graph community games. PLOS ONE , 18(4):1 35. Benati, S., Puerto, J., and Temprano, F. (2024). Modularity for hypergraph clustering: Methodologies and applications. Benders, J. (1962). Partitioning procedures for solving mixed-variables programming problems. Numerische Mathematik , 4:238252. Cheeger, J. (1970). A lower bound for the smallest eigenvalue of the laplacian. Problems in analysis , 625(195-199):110. Cook, S. A. (1971). The complexity of theorem-proving procedures. In Proceedings of the Third Annual ACM Symposium on Theory of Computing , page 151158. Association for Computing Machinery. Dantzig, G. B., Orden, A., and Wolfe, P. S. (1954). Notes on Linear Programming: Part I: The Generalized Simplex Method for Minimizing a Linear Form Under Linear Inequality Restraints . RAND Corporation. Donath, W. E. and Homan, A. J. (1973). Lower bounds for the partitioning of graphs. IBM Journal of Research and Development , 17(5):420425. Everett, H. (1963). Generalized lagrange multiplier method for solving problems of optimum allocation of resources. Operations Research , 11(3):399417. 33
34 Bibliography Fiedler, M. (1975). A property of eigenvectors of nonnegative symmetric matrices and its application to graph theory. Czechoslovak mathematical journal , 25(4):619633. Fischetti, M., Salazar González, J. J., and Toth, P. (2006). The Generalized Traveling Salesman and Orienteering Problems , volume 12, pages 609662. Girvan, M. and Newman, M. E. J. (2002). Community structure in social and biological networks. Proceedings of the National Academy of Sciences , 99(12):78217826. Goldschmidt, O. and Hochbaum, D. S. (1994). A polynomial algorithm for the k-cut problem for xed k. Mathematics of Operations Research , 19(1):2437. Grötschel, M. and Wakabayashi, Y. (1990). Facets of the clique partitioning polytope. Mathematical Programming , 47:367387. Grötschel, M. and Wakabayashi, Y. (1989). A cutting plane algorithm for a clustering problem. Math. Program. , 45:5996. Jonnalagadda, A. and Kuppusamy, L. (2018). A cooperative game framework for detecting overlapping communities in social networks. Physica A: Statistical Mechanics and its Applications , 491:498515. Karp, R. (1972). Reducibility among combinatorial problems. volume 40, pages 85103. Karush, W. (1939). Minima of functions of several variables with inequalities as side conditions. Master's thesis, Department of Mathematics, University of Chicago. Khachiyan, L. G. (1980). Polynomial algorithms in linear programming. USSR Computational Mathematics and Mathematical Physics , 20(1):5372. Kuhn, H. W. and Tucker, A. W. (1951). Nonlinear programming. In Proceedings of the Second Berkeley Symposium on Mathematical Statistics and Probability, 1950 , pages 481492. University of California Press. Lancichinetti, A., Fortunato, S., and Radicchi, F. (2008). Benchmark graphs for testing community detection algorithms. Physical Review EStatistical, Nonlinear, and Soft Matter Physics , 78(4):046110. Markowitz, H. (1952). Portfolio Selection. Journal of Finance , 7(1):7791. Nesterov, Y. and Nemirovskii, A. (1994). Interior-Point Polynomial Algorithms in Convex Programming . Society for Industrial and Applied Mathematics. Nesterov, Yu. and Nemirovski, A. (1992). Conic formulation of a convex programming problem and duality. Optimization Methods and Software , 1(2):95115. Newman, M. (2010). Networks: An Introduction . Oxford University Press, Inc., USA. Newman, M. and Girvan, M. (2004). Finding and evaluating community structure in networks. Physical review. E, Statistical, nonlinear, and soft matter physics , 69:026113.
Bibliography 35 Newman, M. E. J. (2004). Analysis of weighted networks. Phys. Rev. E , (70 (5), 056131). Nickel, S. and Puerto, J. (2009). Location theory. A unied approach , volume 66. Owen, G. (1982). Game Theory . Academic Press. Ponce, D., Puerto, J., and Temprano, F. (2024). Mixed-integer linear programming formulations and column generation algorithms for the minimum normalized cuts problem on networks. European Journal of Operational Research , 316(2):519538. Pothen, A., Simon, H. D., and Liou, K.-P. (1990). Partitioning sparse matrices with eigenvectors of graphs. SIAM journal on matrix analysis and applications , 11(3):430 452. Saldanha-da Gama, F. and Wang, S. (2024). Facility Location Under Uncertainty: Models, Algorithms and Applications . Shapley, L. S. (1951). Notes on the n-person game ii: The value of an n-person game. Shi, J. and Malik, J. (2000). Normalized cuts and image segmentation. IEEE Transactions on Pattern Analysis and Machine Intelligence , 22(8):888905. Sibson, R. (1973). Slink: An optimally ecient algorithm for the single-link cluster method. Comput. J. , 16:3034. Stevens, G. C. (1989). Integrating the supply chain. International Journal of Physical Distribution & Logistics Management , 19:38. Turing, A. (1936). On computable numbers, with an application to the entscheidungsproblem. Proceedings of the London Mathematical Society , 42(1):230265. Wasserman, S. and Faust, K. (1994). Social Network Analysis: Methods and Applications . Structural Analysis in the Social Sciences. Cambridge University Press. Wertheimer, M. (1938). Laws of organization in perceptual forms. Wolfe, P. (1961). A duality theorem for non-linear programming. Quarterly of Applied Mathematics , 19(3):239244. Zhang, S., Wang, R.-S., and Zhang, X.-S. (2007). Identication of overlapping community structure in complex networks using fuzzy c-means clustering. Physica A: Statistical Mechanics and its Applications , 374(1):483490.
Chapter 6 A mathematical programming approach to overlapping community detection 37
Physica A 602 (2022) 127628 Contents lists available at ScienceDirect Physica A journal homepage: www.elsevier.com/locate/physa A mathematical programming approach to overlapping community detection Stefano Benati a,1, Justo Puerto b,1, Antonio M. Rodríguez-Chía c,1, Francisco Temprano b,∗,1 aDipartimento di Sociologia e Ricerca Sociale, Università di Trento, Via Verdi 26, 38122 Trento, Italy bIMUS, Universidad de Sevilla, Avda. Reina Mercedes s/n, 41012 Sevilla, Spain cFaculty of Sciences, Universidad de Cádiz, Avda. República Saharaui, 11510 Puerto Real (Cádiz), Spain article info Article history: Received 25 November 2021 Received in revised form 24 February 2022 Available online 2 June 2022 Keywords: Complex networks Overlapping community detection Modularity Fuzzy membership Mathematical programming abstract We propose a new optimization model to detect overlapping communities in networks. The model elaborates suggestions contained in Zhang et al. (2007), in which overlapping communities were identified through the use of a fuzzy membership function, calculated as the outcome of a mathematical programming problem. In our approach, we retain the idea of using both mathematical programming and fuzzy membership to detect overlapping communities, but we replace the fuzzy objective function proposed there with another one, based on the Newman and Girvan’s definition of modularity. Next, we formulate a new mixed-integer linear programming model to calculate optimal overlapping communities. After some computational tests, we provide some evidence that our new proposal can fix some biases of the previous model, that is, its tendency of calculating communities composed of almost all nodes. Conversely, our new model can reveal other structural properties, such as nodes or communities acting as bridges between communities. Finally, as mathematical programming can be used only for moderate size networks due to its computation time, we proposed two heuristic algorithms to solve the largest instances, that compare favourably to other methodologies. ©2022 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/). 1. Introduction The community detection problem is one of the most studied and interesting problems in networks science. It consists in classifying the units of a population into groups using only information about their links, so that units of the same group can be interpreted as communities. The model is formulated on a network G=(V,E) in which the vertices Vstand for the units, such as individuals, companies, and so on, and the edges Estand for the relations between units, such as kinship, commercial alliances, and so on. Community detection applications can be found in several disciplines, such as biology [1], ecology [2,3], economics [4], sociology [5,6], and many more. A crucial assumption of a standard community detection model is that communities form a partition of V, that is, every unit belongs to one and only one community. This assumption is problematic, as there are applications in which units can realistically belong to two or more communities. For example, the sociological literature documented the role of ∗Corresponding author. E-mail addresses: [email protected] (S. Benati), [email protected] (J. Puerto), [email protected] (A.M. Rodríguez-Chía), [email protected] (F. Temprano). 1All the authors contributed evenly to all the aspects of this research. https://doi.org/10.1016/j.physa.2022.127628 0378-4371/©2022 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license (http://creativecommons. org/licenses/by-nc-nd/4.0/).
S. Benati, J. Puerto, A.M. Rodríguez-Chía et al. Physica A 602 (2022) 127628 individuals who close gaps between communities as members of two communities, see [7], and playing an important role on the functioning of the network as they are bridges, for example, they foster information spreading. Another example is the case of the protein network described in [1], in which proteins are nodes and communities are proteins that carry on one task, but one protein can interact with other communities to accomplish other functions and so it belongs to two or more communities. A standard community detection model may fail to recognize these units, and so emerged the quest for determining what is the best way for finding overlapping communities, as documented in the seminal paper [8] and the survey [9]. The principles that were applied by the algorithms detecting overlapping communities are the same principles used to detect non-overlapping communities, for example constructive methods, Fellows et al. [10], hierarchical clustering, Lancichinetti et al. [11], optimization methods such as mathematical programming in [8], and so on: Other approaches can be found in [12–16] and the survey by Xie et al. [9]. Here, we consider the methods that are based on the comparison of an objective function: That is, given two possible partitions, the best is the one with the highest objective value. When an objective function is used, then there is a strong consensus that the modularity, as defined by Newman and Girvan [17], is the most appropriate measure to detect disjoint communities. Nevertheless, the question of measuring the goodness (or quality) of overlapping communities is more controversial, as different measures were proposed by the literature. In [8], Newman’s modularity function is extended using the fuzzy-cmean, so that a fuzzy membership function accounts for the possibility of multiple communities membership. As that objective function is non-linear, optimal solutions were not calculated and applications are tested using an heuristic algorithm. The contribution in [18] follows the same stream, as a new modularity function is introduced to account for multiple memberships. In this case, group memberships are represented as probabilities instead of fuzzy measures and a genetic heuristic algorithm is proposed to calculate the communities. In [1], the standard Newman’s modularity function is optimized twice: Firstly, to find a node partition in which some nodes could be identified as bridges. Secondly, modularity is applied as an optimization model in which some node can be assigned to more than one community, through an appropriate modification of the mathematical programming model. In [19], a cooperative game is defined on the network and its characteristic function is used to measure players’ Shapley values. Communities are identified as stable coalitions, e.g. the ones in which no player can unilaterally improve its Shapley values by moving into another coalition. Again, no attempt is made to formalize the optimization problem and communities are calculated through heuristics. Unfortunately, the use of heuristic methods to solve a well-posed mathematical programming problem can bring about biases in computing the correct overlapping communities. Indeed, heuristics are devised combining mathematical programming concepts with constructive rule-ofthumb considerations, to the point that empiric solutions can be very far from the optimal. As a matter of fact, when we tested some of these models with more accuracy, we found that sometimes they calculate inconsistent communities, for example, this is the case of the model proposed in [8]. As we documented in our contribution, when we calculated the optimal solution of that mathematical programming problem, we discovered that some overlapping communities can be the same community counted twice. Here, we propose an amendment of the model of Zhang et al. [8], that we could prove it calculates meaningful overlapping communities, e.g., they reveal structural properties of nodes, or of group of nodes, that were not detected by the previous contribution. The new model improves basic ideas contained in [8] with new contributions. They are: •Using a fuzzy membership function uik to determine whether a node ibelongs to a community Ck; •Calculating uthrough the optimization of an objective function; •Using as objective function a variation of the Newman and Girvan’s modularity index. To begin with, our tests revealed that the fuzzy modularity function optimized in [8] biases optimal solutions to grand communities, a.g. communities that are formed by all or almost all nodes. Therefore, we introduce a different fuzzy modular objective function, that avoids this bias. Next, we formulate a Mixed-Integer Linear Programming (MILP) model that maximizes the fuzzy modular function under some linear and integer constraints. The advantage of this approach is that, at least for instances of moderate size, the optimal solution can be calculated exactly with off-the-shelf solvers. Calculating optimal instead of empiric overlapping communities has several advantages, the most important is that the quality and reliability of the communities can be established without the biases that are due to the use of an heuristic. When we applied our new model to some standard benchmark networks, we found meaningful overlapping communities. Unfortunately, MILP problems can be solved for moderate size instances only. However, both the new objective function and the MILP structure lead naturally to new heuristic procedures, that takes advantage of the mathematical formulation of the problem. When applied to large datasets, the heuristics are fast and reliable. After this introduction, the paper is organized as follows. In Section 3we provide an exact formulation of Zhang et al.’s model as Mixed Integer Concave Problem (MICP) so we can calculate its optimal solution for some test problems exactly. We discover that optimal solutions are quite different from the heuristic ones reported in [8] and, unfortunately, they are to a large extent meaningless, as they are the same community counted twice, or they are communities formed by all nodes but one. The misbehaviour is due to the fuzzy modularity index that was initially proposed, so we suggest a way to correct the index to avoid those inconsistent results. Using the new index, a new MILP model is formulated and successfully tested in Section 4. Next, in Section 5, we propose new heuristic algorithms to calculate overlapping communities for large size networks. It can be seen that they find optimal communities in short computing times. The paper concludes with some remarks and suggestions for future research in Section 6. 2
S. Benati, J. Puerto, A.M. Rodríguez-Chía et al. Physica A 602 (2022) 127628 2. The modularity function for overlapping communities Let G=(V,E) be a non-oriented and non-weighted network (or graph), with vertex set V= {1,...,n}and edge set E, represented by the adjacency matrix Aij, e.g., aij =1 if (i,j)∈E,aij =0 otherwise. Let m= |E|and kithe adjacency degree of node i∈V. In [17], modularity optimization is proposed to detect the non-overlapping communities of a graph. The modularity function compares the edge density between nodes of the same community with the expected edge density between the same nodes, but under the assumption that they do not form a community. The expected edge density is obtained from what is called the configuration model. The configuration model is the random graph obtained when edges are placed between nodes randomly, but keeping constant the adjacency degree of each node. In the configuration model, for two nodes iand jthe expected number of edges between them is kikj 2m. Let P= {C1,...,Cq}be a partition of V. Then, the modularity function of the partition P= {C1,...,Cq}is: 1 2m∑ i,j∈V(Aij −kikj 2m)δ(i,j),(1) where δis the Kronecker function, attaining the value 1 if iand jbelong to the same community, 0 otherwise. The formula of the modularity can be extended to weighted graphs as well. The entries Aij are replaced by weights Wij,kiby the sum of the weights associated with arcs adjacent to iand mis replaced by the total sum of weights W=∑(i,j)∈EWij, as described in [20]. Let ncbe the optimal number of communities and Vkbe the set of nodes that belong to community k. The previous modularity function (1) can be expressed as: nc ∑ k=1⎛ ⎝∑i,j∈VkAij 2m−(∑i∈Vkki 2m)2⎞ ⎠.(2) Maximizing the modularity reveals the network community structure, defined as a partition of V. Nevertheless, there are some applications in which a hard partition (hard in the sense that a node can belong to only one community) cannot reveal interactive effects between nodes of different communities. When communities overlap, functions (1) or (2) are not sufficient to determine the network structure, so that, in [8], it is proposed to combine modularity with soft partitions. Instead of assuming that the membership of unit ito community kis represented by a 0-1 number, there is a fuzzy value 0 ≤uik ≤1 that represents the membership of node ito community k. This value uik is called membership function: If uik =1, then ibelongs to community kfor certain, if uik =0, then idoes not belong to community kfor sure, while values uik in the range represent the uncertainty of the membership. It follows that the membership sum is one: ∑nc k=1uik =1∀i∈V, where ncis the number of communities. Membership functions identify the communities structure. For a given threshold value λ, a community Vkis defined as the set of nodes whose membership exceeds the threshold λ,Vk= {i∈V:uik > λ}. From the definition of uand λ, communities Vk,k=1,...,ncmay not form a partition, but are admitted to overlap: A node imay belong to more than one community. Clearly, overlapping communities depend on λ: If λ > 0.5, than nodes belong to one community at most and no community can overlap. Next, if 0.334 < λ ≤0.5, each node can belong to a maximum of 2 communities. In general, assuming pis a positive integer, if 1 p+1< λ ≤1 p, each node can belong to a maximum number pof communities. The fuzzy modularity function for overlapping communities, introduced in [8], is: nc ∑ k=1⎛ ⎝∑i,j∈VkAij uik+ujk 2 2m−(∑i,j∈VkAij uik+ujk 2+∑i∈Vk,j/∈VkAij uik+(1−ujk) 2 2m)2⎞ ⎠.(3) This function is a modification of the original modularity function (2). It is obtained by weighting each edge (i,j) by the average of uik and ujk if i,j∈Vk, or by the average of uik and 1 −ujk if i∈Vk,j∈ Vk. Note that the expression (3) can be applied to weighted graphs, too. It is sufficient to change the entry Aij with an edge weight, Wij, and mreplaced by the total sum of weights W=∑(i,j)∈EWij. 3. An exact mathematical programming formulation of the maximum fuzzy modularity problem Finding the ncoverlapping communities that maximizes the function (3) is a hard problem to solve. Indeed, it is a mixed binary polynomial optimization problem which can be easily proven to be NP-hard [21] and, in [8], only a heuristic procedure is proposed indeed. Unfortunately, heuristic procedures may result with sub-optimal solutions that can be very different from the optimal. Therefore, to test the effectiveness of that model, we formulated a mathematical programming model that exactly maximizes function (3). Problem variables are the membership functions uik, such that 0 ≤uik ≤1 for all iand k. Next, hard membership functions xik, for all iand k, that depend on variables uik and threshold parameter λ, are defined in the following way: xik ={1,if node iis assigned to community k, that is, if uik > λ, 0,otherwise. 3
S. Benati, J. Puerto, A.M. Rodríguez-Chía et al. Physica A 602 (2022) 127628 Fig. 10. Zachary’s Karate club communities maximizing Problem (NEW-MOD) with nc=4 and λ=0.25. Fig. 11. Highland tribes communities maximizing the Newman and Girvan’s modularity. To summarize the previous tests, there is some flaw in the definition of optimal overlapping communities, as proposed in [8], that regards the formulation of the objective function (3). The previous contribution could not recognize the inconsistency, because optimization was implemented through the use of an heuristic procedure. The heuristic procedure stacked in favour of sub-optimal solutions that were actually so far from optimal ones that they confuse reasonable communities with optimal. To motivate this argument, note that a set of overlapping communities correspond to fixing hard assignment variables xto 0 or 1. Then, with xbeing fixed, problems (F-MOD) and (F-MOD-NI) can be solved to calculate the corresponding soft assignment variables u. Finally, the objective function (3) can be calculated and compared using various x. Table 1reports, the objective function for: (i) overlapping communities calculated by problem (F-MOD), (ii) overlapping communities calculated by problem (F-MOD-NI), (iii) overlapping communities calculated by the heuristic procedure in [8] and reported in Figs. 2 and 3. As can be seen, the objective function of the suboptimal solutions are very far from the optimal ones. 10
S. Benati, J. Puerto, A.M. Rodríguez-Chía et al. Physica A 602 (2022) 127628 Fig. 12. Highland tribes communities maximizing Problem (NEW-MOD) with nc=3 and λ=0.25. Fig. 13. Optimal Highland tribes communities maximizing (NEW-MOD) with nc=4 and λ=0.4. 11
S. Benati, J. Puerto, A.M. Rodríguez-Chía et al. Physica A 602 (2022) 127628 Fig. 14. Zebra communication community structures. 4. A new fuzzy modularity function for overlapping community detection Our previous experiments showed that formula (3) is not appropriate to detect overlapping communities, as it leads to optimal solutions composed of almost all the vertices. The reason to this is that, if node i∈Vkand j/∈Vk, weighting the edges between iand jby the average of uik and 1 −ujk introduces a bias in the sum. Indeed, if j∈ Vk, then (1 −ujk) is large, but this term appears as a subtraction in (3), and therefore this is an incentive to create large communities, just to avoid these negative terms. Therefore, a reasonable transformation of formula (3) could be to exclude those negative terms and to retain only the positive ones, [24,25]. To formulate the new measure, we introduce some hypothesis, as those that inspired the modularity function in [17]. There, similarity between units are compared to the ones that were obtained randomly, e.g., the ones of the so-called configuration graph. In our proposal, instead of generating one random graph, we generate ncrandom graphs and we assume edge weights calculated as the average of their membership for each of the nccommunities. Therefore, the new objective function is: nc ∑ k=1∑ i,j∈Vk(Aij −kikj 2m)uik +ujk 2,(30) 12
S. Benati, J. Puerto, A.M. Rodríguez-Chía et al. Physica A 602 (2022) 127628 Fig. 15. Windsurfers communities maximizing Newman and Girvan’s modularity are the same maximizing Problem (NEW-MOD) with nc=2 and λ=0.4,0.25. where the condition i∈Vkcorresponds to the hard assignment resulting from uik ≥λ. Note that the objective function (30) is linear, while the objective function (3) is quadratic, leading to a mixed-integer optimization problem that should be solved faster than its quadratic counterpart. Moreover, Eq. (30) can be extended to consider weighted graphs too. Entry i,jof the adjacency matrix, Aij is replaced by Wij,ki, and the weights sum is W=∑(i,j)∈EWij. The new model is: (NEW-MOD) max 1 2m nc ∑ k=1 n ∑ i,j=1(Aij −kikj 2m)wijk +wjik 2(31) s.t.: (4),(5),from (8) to (13),(20),(22),(23). In Problem (NEW-MOD), two parameters appear as input: They are the membership threshold λand the total number of communities nc. There is not an optimal choice for these parameters, rather, there is a trade-off between them and the value of the objective function, as it happens similarly for the k-means model. The correct choice depends on the application at hand. Here, we propose to solve maximum modularity with disjoint community to calculate nc, and then using it as the input of Problem (NEW-MOD) to let these communities overlap, e.g. controlling whether the possibility of multiple memberships increases the objective function. Nevertheless, this choice is not exclusive and other rule-of-thumb can be used to the purpose. Next, parameter λcontrols for the maximum number of communities to which a node can belong, for example and as discussed previously, for λ∈ [0.34,0.50) a node can belong to two communities at most. Therefore, admissible values are λ∈ {1 2,1 3,..., 1 nc}. Nevertheless, choosing the right one depends on the application. After having implemented formulation (NEW-MOD) in Python and solved with the Gurobi’s solver, we applied it to the first example in [8] and on the Zachary’s karate club network. Our results are reported in Figs. 5 and 6(a), respectively. As can be seen, for the first example, Fig. 5 replicates the communities that have been found in [8]. For what concerns the karate club, we can compare Fig. 2 with 6(a) and observe that the results of the new model are communities similar to the ones that were calculated with a constructive algorithm in [8]. To see the effect of varying λ, the results with λ=0.1 are reported in Fig. 6(b). Similar overlapping communities can be seen, but more intersection nodes appear due to the smaller value of λ. One important advantage of using mathematical programming to formulate and solve the overlapping community detection is that additional features or constraints that one expects from communities can be explicitly modelled as linear inequalities of the problem constraints. For example, as suggested in [1], a researcher may know from other qualitative sources that some nodes, e.g. actors of the networks, are acting as bridges between groups. In this case, these nodes should be included in at least two communities from the beginning, that is, from the problem formulation. This property can be imposed over a node iby the inequality: nc ∑ k=1 xik ≥2 in the constraints of Problem (NEW-MOD). Another possibility is to require that full inclusion among overlapping communities is explicitly forbidden. This can be enforced including in the formulation inequalities (25)–(29). 13
S. Benati, J. Puerto, A.M. Rodríguez-Chía et al. Physica A 602 (2022) 127628 Fig. 16. Zachary’s karate club community structure obtained by the Algorithm 1with λ=0.25. Fig. 17. Highland tribes community structure obtained by the Algorithm 1with λ=0.25. 14
S. Benati, J. Puerto, A.M. Rodríguez-Chía et al. Physica A 602 (2022) 127628 Fig. 18. Zebra communication community structure obtained by the Algorithm 1. Table 1 Fuzzy (3) comparisons of the most outstanding structures. Dataset ncλF-MOD F-MOD-NI [8] Zachary’s karate club 3 0.25 0.667 0.65 0.445 American college football team 10 0.1 0.9 0.887 0.619 As an example of the previous remark, we calculated the structure of the Zachary’s karate club network with nc=3 and λ=0.25 but imposing some nodes to belong to two communities, acting as bridges. Following the solution reported in [8], we let nodes 1,9,10 and 31 belong to two communities. The solution obtained is depicted in Fig. 7. The same 15
S. Benati, J. Puerto, A.M. Rodríguez-Chía et al. Physica A 602 (2022) 127628 Fig. 19. American college football team community structure obtained by the Algorithm 1with λ=0.1. overlapping structure as in Fig. 6(a) are obtained, but now nodes 9,10 and 31 are at the intersection of the communities black and grey. Applied to the American football instance, after 24 h of computing time, the linear solver could not certify the optimality of its incumbent solution to problem (NEW-MOD). This incumbent solution is depicted in Fig. 8. Even though the solution has not been proved to be the best, still its objective function, (3) e.g f=0.632, is higher than the objective function of the solution calculated in [8]. Moreover, it can be observed in Fig. 8 that there are many intersection nodes, as λhas a low value, communities tend to overlap more. In the following, we apply formulation (NEW-MOD) using the following procedure: First, we calculate the optimal disjoint communities that maximizes the Newman and Girvan’s modularity to determine the parameter nc, the number of these communities. Then we apply Problem (NEW-MOD) with varying ncand λas input parameters. The first example is again Zachary’s karate club. The optimal non-overlapping communities are reported in Fig. 9, the output is nc=4. When admitting overlapping communities, the results are reported in Fig. 10. As can be seen, overlapping communities identify nodes that are bridges between communities, as their edges are adjacent to different groups, as is the case of nodes 1, 3, 12, 24, 34. This result shows that overlapping communities provide important information about the structural properties of nodes, information that is not available when communities are disjoint. The second example is the alliances network between Highland tribes of New Guinea, Read [26]. The optimal non-overlapping communities are reported in Fig. 11, the output is nc=3. When admitting overlapping communities, and with λ=0.25, the results are reported in Fig. 12, while, for nc=4 and λ=0.4, the results are reported in Fig. 13. The detected communities contain a high internal edge density and the intersection nodes share connections with various communities. These results can be interpreted easily and robust to different parameters combination The third example is the zebra communication network, [27]. The optimal non-overlapping communities are reported in subfigure (a) of Fig. 14, the output is nc=4. 16
S. Benati, J. Puerto, A.M. Rodríguez-Chía et al. Physica A 602 (2022) 127628 Table 2 Computational results of the solution methods. Dataset ncλSolving method Time (s) Objective value Zachary’s karate club 4 0.25 Exact model 42315 0.442831 Algorithm 1.i. 4.26 0.441979 Algorithm 1.ii. 0.66 0.440787 Algorithm 1.iii. 0.4 0.440787 Zachary’s karate club 3 0.25 Exact model 11573 0.41415 Algorithm 1.i. 131 0.41415 Algorithm 1.ii. 14.68 0.41415 Algorithm 1.iii. 5.47 0.41415 Highland tribes 3 0.25 Exact model 23715 0.191439 Algorithm 1.i. 0.43 0.191439 Algorithm 1.ii. 0.17 0.191439 Algorithm 1.iii. 0.08 0.184379 Zebra communication 4 0.4 Exact model 1463 0.282266 Algorithm 1.i. 0.72 0.282266 Algorithm 1.ii. 0.22 0.282266 Algorithm 1.iii. 0.19 0.282266 Zebra communication 4 0.25 Exact model 3682 0.284342 Algorithm 1.i. 1.33 0.284342 Algorithm 1.ii. 0.35 0.284342 Algorithm 1.iii. 0.15 0.282911 American college football team 10 0.1 Exact model 86400 0.619 Algorithm 1.i. 1902 0.6345487 Algorithm 1.ii. 136 0.6345 Algorithm 1.iii. 56.78 0.616872 When admitting overlapping communities, and testing for λ=0.25,and 0.4, it can be seen that as λis smaller, there are more nodes in the intersections between communities. The fourth example is the windsurfers network, [28]. The optimal nonoverlapping communities result in nc=2. When admitting overlapping communities, still disjoint communities are detected, see Fig. 15. This is because there are two distinguished communities and the model cannot find nodes behaving as bridges between groups. From the four applications we can conclude that: •Overlapping communities provide additional information about the structure of the connection between groups, for example, identifying nodes interpreted as bridges. •Varying parameter λis an effective tool to let communities of different shape emerge, identifying what nodes are also the most influential, as they belong to more than two communities. •The overlapping model is flexible enough to guarantee a non-overlapping community as the outcome, when data suggest so. 5. Heuristic algorithms to approximate overlapping communities If λis fixed to 1, then Problem (NEW-MOD) is the maximum modularity problem, a problem that is known NP-hard, therefore NP-hard itself: They are problem for which a polynomial time algorithm does not exist, unless P = NP. In practice, this negative result implies that optimal solutions can be obtained only for instances of moderate size. As can be seen in the application to the American college football teams this size is of the order of a few tens. Nevertheless, the MILP formulation can be used to obtain good approximations of optimal solutions through the use of heuristic. The first algorithm that we propose is based on local search. It consists of an iterative method applied to feasible solutions that are improved by modifying some value, e.g. membership functions or hard assignments, until no improvement is possible. In that case, we say that a local optimum has been reached. More formally, let ncbe the maximum number of communities and λbe the threshold of the membership value, let Π= {V1,...,Vn∗ c}be a set of n∗ ccommunities of V,n∗ c≤nc. We say that Πis a feasible solution if: (1) every node belongs to at least one community, that is, ⋃n∗ c k=1Vk=V, (2) no community is a subset of a larger one, that is, ∄k,r=1,...,nc,k= r, such that Vk⊆Vr, and (3) there is a feasible membership solution ufor the partition Π, i.e. ∀i∈Vthe inequality 1 |{k=1,...,n∗ c:i∈Vk}| ≥λis fulfilled. Assume that Πis a feasible solution. To improve the objective function we consider three types of changes of hard assignments: (1) to add a node to a community or (2) to remove a node from a community, and (3) to swap two nodes 17
S. Benati, J. Puerto, A.M. Rodríguez-Chía et al. Physica A 602 (2022) 127628 Table 3 Computational results comparing Algorithm 2with Clique percolation and Fuzzy c-means. nλNumber of communities Solving method Time (s) Objective value 333 0.5 13 Large scale algorithm 0.84 0.4381 15 Fuzzy c-means 63.35 0.36 172 Clique percolation 26.67 0.3135 333 0.2 13 Large scale algorithm 0.84 0.44 13 Fuzzy c-means 64.32 0.3614 109 Clique percolation 27.98 0.3871 747 0.5 6 Large scale algorithm 7 0.526 9 Fuzzy c-means 262.59 0.506 – Clique percolation – – 747 0.2 6 Large scale algorithm 7.28 0.5286 9 Fuzzy c-means 270.174 0.512 – Clique percolation – – 224 0.5 5 Large scale algorithm 0.51 0.2691 12 Fuzzy c-means 37.33 0.254 100 Clique percolation 1102.26 0.1679 224 0.2 5 Large scale algorithm 0.91 0.2933 12 Fuzzy c-means 38 0.26169 90 Clique percolation 1268.74 0.1893 534 0.5 10 Large scale algorithm 3.9 0.62656 10 Fuzzy c-means 140.74 0.63 246 Clique percolation 70.68 0.488 534 0.2 10 Large scale algorithm 4.5 0.6445 10 Fuzzy c-means 149.83 0.635 93 Clique percolation 77.8 0.6224 1034 0.5 6 Large scale algorithm 39.36 0.5357 7 Fuzzy c-means 499.18 0.52428 – Clique percolation – – 1034 0.2 6 Large scale algorithm 46.9 0.5401 7 Fuzzy c-means 519.15 0.52689 – Clique percolation – – from two different communities. These moves can be applied if and only if the new obtained solution is feasible. For instance, the first and the second movement cannot be applied if it results in some inclusion between communities. More formally, for i∈Vand k=1,...,nc, let the triplet (i,k,1) be the move of adding node ito community kand let the triplet (i,k,2) be the move of removing node ifrom community k, let the 5-tuple (i,k,i′,k′,3) be the move of swapping nodes iand i′between communities kand k′respectively. Assume that Πis a feasible solution, represented by hard assignments xand membership functions u, then, to calculate the objective function (30), membership functions umust be determined too: ucan be calculated in the following ways: (i) Exact calculation of u: use formulation (NEW-MOD) to calculate the optimal objective function (30) and its corresponding u. (ii) Approximate calculation of u: keep fixed the membership functions ucorresponding to unchanged assignments x, and find an approximate value uik only for the assignments xik that were modified using formulation (NEW-MOD). This option is less accurate but also reduces complexity. (iii) Approximate calculation of u: Approximate uas follows. If xik =0 then uik =0, while uik =p, with pa constant term, for the value for which xik =1, so p=1 ∑n∗ c j=1xij . This is the least accurate approximation, but also the simplest and fastest, as it does not require any optimization. Finally, the interchange heuristic is applied to the initial solution xcalculated maximizing modularity (1) with nonoverlapping solution. Further diversification can be obtained by choosing the initial solution xrandomly, e.g. using the so-called random restart. The pseudo-code of the algorithm is summarized in Algorithm 1: 18
S. Benati, J. Puerto, A.M. Rodríguez-Chía et al. Physica A 602 (2022) 127628 Algorithm 1 Heuristic algorithm procedure Local Search for Overlapping Communities Π= {V1,...,Vnc} ← Initial_Subdivision ▷Πfeasible subdivision obtained randomly or by another procedure f←Extended_Modularity(Π)▷The extended modularity (30) is computed by three different ways: i), ii) or iii). local_opt =FALSE ▷Condition for a local optimum while local_opt =FALSE do ∆←Feasible_Moves(Π)▷∆: list of admissible movements for Π. for (i,k,d) in ∆do if d=1 then Vk←Vk∪ {i} δikd ←Extended_Modularity(Π) Vk←Vk\ {i} end if if d=2 then Vk←Vk\ {i} δikd ←Extended_Modularity(Π) Vk←Vk∪ {i} end if end for for (i,k,i′,k′,3) ∈∆do if d=3 then Vk←Vk∪ {i′} \ {i} Vk′←Vk′∪ {i}\{i′} δiki′k′d←Extended_Modularity(Π) Vk←Vk\ {i′}∪{i} Vk′←Vk′\ {i}∪{i′} end if end for (i∗,k∗,d∗)∈arg max{δikd|(i,k,d)∈∆} ▷ Select the move that increases the most (i∗,k∗,i′∗,k′∗,d∗)∈arg max{δiki′k′d|(i,k,i′,k′,d)∈∆} ▷ Select the move that increases the most if δi∗k∗i′∗k′∗d∗>max{f, δi∗k∗d∗}then f←δi∗k∗i′∗k′∗d∗▷Update f Vk∗←Vk∗∪ {i′∗}\{i∗} Vk′∗ ←Vk′∗ ∪ {i∗}\{i′∗} ▷ Update Π else if δi∗k∗d∗>fthen f←δi∗k∗d∗▷Update f if d∗=1then Vk∗←Vk∗∪ {i∗} ▷ Update Π else Vk∗←Vk∗\ {i∗} ▷ Update Π end if else local_opt =TRUE end if end if end while return Π▷Return the local optimum end procedure If the initial solution Πis obtained randomly, then the interchange can be repeated for a maximum of tmax initial solution. For each attempt t, we obtain local optimal objective function ftand overlapping communities Πt. Finally, approximate solution of problem (NEW-MOD) is the best local optimum. In test problems, we run the algorithm with the three different methods to calculate membership function u. To assess their quality, we applied Algorithm 1to the initial solution Πcalculated by the maximum modularity (1), as in this way we prevent potential biases caused by random starting solutions. The communities found in this way are reported in Figs. 16,17,18 and 19. As can be seen, when the membership functions uare calculated exactly, e.g. by solving the optimization problem, then the heuristic communities are the optimal. While, when uare approximated, the heuristic communities only differ for very few nodes to the optimal ones. Next, in Table 2, we compared Algorithm 1with the three variants for computing uwith the optimal solution of problem (NEW-MOD). It can be seen that the heuristic algorithms reduce the computational time at the cost of decreasing the objective function only to a small amount. When the optimal solution is not available, such as the case of the American football data, two of the heuristic algorithms could obtain a better objective value than the MILP truncation after 24 h of computing. For the Zachary’s karate club case with nc=3, since the optimal number of communities in the non-overlapping case is nc=4 and λ=0.25, we do not use this as an initial solution. Instead, we perform a multistart strategy with ten iterations starting with 3 randomly chosen communities. Algorithm 1works well when input data are small or medium sized networks, e.g. networks with hundreds of nodes, but computational times might be too high for large sized networks, for example instances with more than 1000 nodes. For this reason, we developed a variation of the previous Algorithm, reported in Algorithm 2, in which some operations 19
The objective function is merely a simple statistic that evaluates partitions or node assignments. As such, it can be used to compare alternative community structures and to decide what is the most meaningful. One of the most popular statistic is modularity, see [1]. Modularity is an index that, for a given partition, compares the arc density of a subset with the one that is obtained on the assumption of node random pairings. The highest the modularity, the most connected are the nodes within a community, allowing a clear substantial definition of what is a community. The extension of the modularity to the case of overlapping communities has been proposed in [9], using fuzzy membership functions that are optimized using the fuzzy-cmeans algorithm. This method has been elaborated further in [10–13], where the standard modularity function is modified by node or arc weights, representing node affinity, fuzzy memberships, or other. An alternative version of the objective function proposed in [9] is presented in [14], fixing some biases of the original one. In [7], it is proposed to maximize the modularity function, but with some additional constraints that allow some nodes to belong to more than one community. These nodes are referred to as bridges. In [15], communities are defined as stable coalitions of a cooperative game. In a cooperative game, a coalition is stable if every member does not take any advantage in leaving the coalition to obtain a better payoff elsewhere, so a community is based on the concept of a common interest. There is a large room to define this common interest through any game characteristic function, such as market, voting, matching games, and so on. To just consider the topological network properties, such as the arc density and the node common neighbors, in [15] a weighted graph community game is proposed, with arc weights defined on some peculiar topological indicators. Next, an objective function is proposed to discern between alternative community structures and a constructive heuristic is implemented to find them. In our contribution, we formulate the problem of finding communities as stable coalitions proposed in [15], as a mixed-integer linear programming problem. In this way, taking advantage of existing software, we can calculate the optimal communities of that model without resorting to any heuristic consideration. As a result, we can evaluate the optimal solutions of that model without the biases due to the use of the heuristic. Indeed, we found that the communities proposed in [15] are far from the optimal ones and, unfortunately, optimal ones are inconsistent too, in the sense that they do not correspond to what empirically one expects to find out. As it will be discussed, we argue that the reason of the inconsistency is on how costs of the weighted graph community game are defined and therefore we proposed a correction to them. Our correction follows the spirit of the modularity function, [1], in which an actual value of a statistic is compared to an expected value in absence of any community structure. We will show that our correction is reliable and effective as, after many computational tests, we showed that our method can recognize the hidden community structure of the networks. As a by-product of our contribution, we note that our cost definition relies on the calculation of the expected value of some network statistics on the assumption that no community is embedded in the network. To have an accurate cost estimate, we elaborated a new theorem to calculate the exact value of these statistics and it is worth to note that this theorem may have an autonomous interest for other applications in which some exact probabilities can be applied, as the same seminal paper [1]. To summarize, the contributions of our paper are the following: 1. We provide a mathematical formulation of the method proposed by [15] to detect the overlapping communities of a network. 2. We show that the communities obtained with this methodology are not the real communities embedded in the network, but we proposed an amendment to the game cost function that correct the bias. PLOS ONE Overlapping communities detection through weighted graph community games PLOS ONE | https://doi.org/10.1371/journal.pone.0283857 April 4, 2023 2 / 35 Puerto and Francisco Temprano acknowledge financial support by the Spanish Ministerio de Ciencia y Tecnologı ´a, Agencia Estatal de Investigacio ´n and Fondos Europeos de Desarrollo Regional (FEDER) with grant number: PID2020114594GB-C21, and Junta de Andalucı ´a with grant number: P18-FR-1422. The authors Stefano Benati, Antonio Manuel Rodrı ´guez-Chı ´a and Justo Puerto also acknowledge partial support from: NetmeetData: Ayudas Fundacio ´n BBVA a equipos de investigacio ´n cientı ´fica 2019 with reference "COMPLEX NETWORK". The author Antonio Manuel Rodrı ´guez-Chı ´a also acknowledges the European Regional Development Fund via projects with grant numbers: FEDER-UCA18-106895 and TED2021-130875B-I00; and the Spanish Ministerio de Ciencia y Tecnologı´a, Agencia Estatal de Investigacio ´n with grant number: PID2020114594GB-C22. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. Competing interests: The authors have declared that no competing interests exist.
3. We propose a heuristic algorithm that can calculate the optimal communities when the exact method fails because of the network size. 4. We apply our new mathematical model to real and artificial test problems and we show its effectiveness and reliability. The paper is organized in 4 sections. In the Introduction, we motivate the paper purpose and summarize its contribution. In Material and methods Section, we formally introduce the overlapping community detection problem and the methods proposed by [15]. There, we design the exact optimization model and observe the finding of inconsistent communities. In Subsection called Detecting overlapping communities as stable coalitions of a cooperative game, we propose an alternative definition of the costs of the weighted graph community game that leads to a different objective function of the optimization model. In Local Stability Exploration Subsection, we present a heuristic algorithm for solving our model for the cases in which the network size is too large to compute the exact solution in a reasonable amount of time. In Results and discussion Section, we compare the exact and heuristic algorithm and then we report some computational results of a controlled experiment on graphs generated according the method proposed in [16] and we show that our method recovers correctly the community structure. The paper ends with some concluding remarks and outlines for future research in the final section, namely Conclusion. Material and methods Detecting overlapping communities as stable coalitions of a cooperative game In [15], a cooperative game on a weighted graph is defined to characterize overlapping communities. The nodes of a graph are considered as the players of a network game, and then the Shapley value is used to characterize stable coalitions, e.g. subsets of nodes in which no player has any incentive to leave. Specifically, the cooperative game (V,φ) is defined on the weighted graph G= (V,E), with V= {1, . . .,n}, e.g. players are nodes labeled from 1 to n, weights W ij (�0) are defined for any edge (i,j)2E, then the game characteristic function is: φðSÞ ¼ X i;j2S i<j Wij;for S�V:ð1Þ That is, the value of coalition Sis the weights sum of the edges of the subgraph induced by S. The model has been called Weighted Graph Community (WGC) Game in the aforementioned paper. When a coalition S�Vis going to form, then the members i2Scan calculate the gain that they can get from it, e.g. what is their share of the payoff φ(S) that they can receive. A standard result of cooperative games is that the share that they can get is the Shapley value of the game restricted to S: For player iand coalition S,i2S, the Shapley value is: φiðSÞ ¼ 1 2X j2S j6¼ i Wij: Hence, the profit of player ifrom coalition Sdepends on the total weight of its connection with the other members of S. PLOS ONE Overlapping communities detection through weighted graph community games PLOS ONE | https://doi.org/10.1371/journal.pone.0283857 April 4, 2023 3 / 35
In [15], a coalition is defined stable if no member of Stakes advantage from swinging from coalition Sto coalition V\S. In mathematical terms it occurs if and only if: φiðSÞ � φiððVnSÞ[figÞ;8i2S:ð2Þ Actually, there are different definition of stable coalitions that can be found in the literature: Stable coalition structures are defined in [17,18], while in [19,20], condition (2) is called the internal stability property. Moreover, in the latter notion of stability, an additional property is imposed requiring that a coalition Sis stable if no member of Stakes advantage from swinging from Sto any other subset S0contained in V\S. This can be formalized as: φiðSÞ � φiðS0[figÞ;8i2S;8S0�VnS:ð3Þ However, we are not developing this issue further and we will remain with definition (2). Formulating a WGC game allows a formal definition of what are the feasible overlapping communities of a network: As a node can belong to more than one stable coalition, communities can overlap. However, a crucial feature of the model is the way in which weights W ij are defined. In [15], the following formula is proposed: Let k i be the adjacency degree of node i (e.g. the number of nodes to which iis connected through an arc), let Pij ¼1 kiþ1 kjbe defined as the partition ratio and let CN ij = (|common neighbors of i and j| + 1)P ij be defined as the neighbourhood ratio of i,j2V, then the weight of the arc (i,j), i6¼jis Wij ¼ CNijPij 4;if ki�1;kj�1and ði;jÞ=2E; Pij;if ki¼1or kj¼1and;ði;jÞ 2 E; 2CNij þPij;if ki>1and kj>1;and ði;jÞ 2 E; 0;otherwise: 8 > > > > > > > < > > > > > > > : ð4Þ The formula was proposed in [15] to consider the node similarity as dependent on both the direct and indirect links between iand j. It is straightforward to observe that W ij �0, but this property has important consequences on the structure of the stable coalitions, as it will be discussed later. For the moment, we focus in the methodology to find all the stable coalitions of a networks. While in [15] a constructive method is proposed, that is, an heuristic technique with some ad-hoc adjustment to find stable coalitions, here we propose a mathematical programming approach in which all considerations about stability discussed in [15] are translated into an objective function and mathematical constraints. We will show that stable coalitions can be represented by linear constraints involving binary variables and then, using an appropriate objective function, stable coalitions can be determined by linear programming. Let n c be the maximum number of communities to which a node can belong to (this is not a binding constraint to the model, since n c can be large enough to include all the feasible stable communities). For i= 1, . . .,nand k= 1, . . .,n c , the model variables are: xik ¼ 1;if node ibelongs to community=coalition Sk; 0;otherwise: 8 > > > < > > > : PLOS ONE Overlapping communities detection through weighted graph community games PLOS ONE | https://doi.org/10.1371/journal.pone.0283857 April 4, 2023 4 / 35
For any i,j= 1, . . .,nsuch that i<jand k= 1, . . .,n c : zijk ¼ 1;if nodes iand jboth belongs to community=coalition Sk; 0;otherwise: 8 > > > < > > > : The relationship between xand z-variables is given by the logical/quadratic constraints z ijk =x ik x jk for all i,j2V,i<jand all k= 1, . . .,n c . Then, the quadratic constraint can be replaced by the linear constraints: zijk �xik;8i;j¼1;. . . ;n;i<j;k¼1;. . . ;nc;ð5Þ zijk �xjk;8i;j¼1;...;n;i<j;k¼1;. . . ;nc;ð6Þ xik þxjk zijk �1;8i;j¼1;...;n;i<j:k¼1;. . . ;nc:ð7Þ Next, using binary x-variables, the stability condition (2) can be characterized by linear constraints too. First, for fixed iand k, consider the quadratic inequality: xik X n j¼1 j6¼ i xjkWij X n j¼1 j6¼ i 1xjk � �Wij 0 B B @1 C C A�0: If x ik = 1, then ibelongs to coalition S k , so that S k must be stable. For the stability, i-player’s Shapley value from coalition S k must be greater than its Shapley value from the opposite coalition (V\S k )[{i}. The term Pj¼1 j6¼ i xjkWij is the Shapley value of coalition S k , as all j’s such that x jk = 1 are all the other players of coalition S k . Conversely, all other j’s such that (1 −x jk ) = 1 are the players excluded from S k . Consequently, Pn j¼1ð1xjkÞWij is the Shapley value of the opposite coalition, (V\S k )[{i}. Finally, their difference must be greater than or equal to 0 for S k to be stable. Next, the above quadratic inequality can be simplified to the following linear one: X n j¼1 j6¼ i xjkWij �Pn j¼1 j6¼ i Wijxik 2;8i¼1;. . . ;n;k¼1;...;nc:ð8Þ Next, it must be imposed that overlapping coalitions/communities must have non-empty difference, e.g. the same coalition is not selected more than once (a coalition must not be contained in a different one). To prevent inclusion, additional variables hare introduced for i= 1, . . .,nand pairs k,rsuch that 1 �k<r�n c : hikr ¼ 1;if ibelongs to community Srand not to community Sk; 0;otherwise: 8 > > > < > > > : PLOS ONE Overlapping communities detection through weighted graph community games PLOS ONE | https://doi.org/10.1371/journal.pone.0283857 April 4, 2023 5 / 35
The relation between xand h-variables is given by the quadratic constraint: h ikr =x ir (1 − x ik ), that can be replaced by three linear constraints as done for z-variables in expressions (5)– (7). To prevent the inclusion of S r in S k , it must be that: X n j¼1 hjkr �xir;81�k<r�nc;8i¼1;. . . ;n:ð9Þ The constraint is binding when x ir = 1. In that case, coalition S r must contain at least one element jthat is contained in S r but not in S k , guaranteeing that S r ⊄S k . To conclude, we introduce inequalities to avoid symmetrical solutions too. Symmetric solutions decrease the efficiency of the Integer Linear Programming solver, as the same structural solution can be obtained by multiple assignments to variables x,z,h, simply giving different labels to coalitions. Note that constraints (9) avoid to replicate the same coalition, so that it is sufficient that, after ranking the communities from the largest to the smallest, they are assigned to decreasing labels k. The following constraints do the task: X n i¼1 xik �X n i¼1 xi;kþ1;8k¼1;...;nc1:ð10Þ Every stable coalition corresponds to a point of the polytope described by the equations and inequalities described so far. To determine what are the most meaningful overlapping communities, in the objective function it is used the nodes Shapley value. If a coalition S k is established, then player i’s Shapley value from coalition S k is: Pn j¼1zijkWij. Therefore, for a set of overlapping communities S k ,k= 1, . . .,n c , the total Shapley value of a player iis the sum of the values it gets from every coalition, that is: X nc k¼1X n j¼1 zijkWij:ð11Þ In [15], the most important overlapping coalitions are determined by maximizing the sum of the Shapley values of all nodes. Therefore, this index will be used as the objective function of the following integer programming formulation: ðFShJKÞmax X nc k¼1X n1 i¼1X n j¼iþ1 zijkWij ð12Þ s.t.: (5)–(10), X nc k¼1 xik �1;8i¼1;. . . ;n;ð13Þ hikr �1xik;8i¼1;...;n;k;r¼1;. . . ;nc;k<r;ð14Þ hikr �xir;8i¼1;...;n;k;r¼1;...;nc;k<r;ð15Þ PLOS ONE Overlapping communities detection through weighted graph community games PLOS ONE | https://doi.org/10.1371/journal.pone.0283857 April 4, 2023 6 / 35
xir xik hikr �0;8i¼1;. . . ;n;k;r¼1;. . . ;nc;k<r;ð16Þ xik 2 f0;1g;8i¼1;. . . ;n;k¼1;. . . ;nc;ð17Þ zijk 2 ½0;1�;8i;j¼1;...;n;i6¼ j;k¼1;...;nc;ð18Þ hikr 2 ½0;1�;8i¼1;. . . ;n;k;r¼1;. . . ;nc;k<r:ð19Þ The objective function (12) represents the sum of the Shapley values for all nodes and communities. Constraints (13) guarantee that every node belongs to at least one community. Constraints (14)–(16) are the linear representations of the h-variables. Finally, constraints (17) define binary variables. Note that in (18) and (19), we can relax the z−and h−variables to be continuous, since the constraints on the x-variables force both to be binary. F Sh−JK is the exact Integer Programming formulation of the model proposed in [15]. However, in that seminal paper the overlapping communities were computed through a heuristic constructive procedure, in which the search for optimal solutions is combined with various ad-hoc adjustments to induce sufficient diversification of coalitions. The advantage of Integer Programming is that the output coalitions of F Sh−JK are exactly the optimal ones, without any bias due to constructive rule-of-thumb procedures. As we will see, this allows us to point out a drawback of the game definition and to suggest a method to adjust it. We apply formulation F Sh−JK , to the Zachary’s karate club network, fixing n c = 3. Optimal overlapping communities can be seen in Fig 1. As can be seen, selected communities are the grand coalition (all the nodes belong to the same coalition) except one node. That is, communities are subsets Ssuch as |S|=n−1, in which the discarded node is the one with less connections. It is hard to believe that those sets are of some interest to researchers, as they are far from the communities that were often identified in the Zachary’s network. The same occurs Fig 1. Zachary’s karate club structure obtained by F Sh−Jk with n c = 3. (a) Community 1, (b) Community 2, (c) Community 3. https://doi.org/10.1371/journal.pone.0283857.g001 PLOS ONE Overlapping communities detection through weighted graph community games PLOS ONE | https://doi.org/10.1371/journal.pone.0283857 April 4, 2023 7 / 35
with all the other problems we tested: Overlapping communities are the grand coalition except one node. The reason of this disappointing result is not the solution method, e.g. exact vs heuristic, or the community definition, e.g. using cooperative games and the Shapley value. Rather, the reason is the way in which weights Ware formulated in (4). As recognized in [15], if W ij � 0 for all i,j, then the cooperative game (V,φ) is convex, that is for two coalitions S,Tsuch that S�Tand i=2T, it always occurs that: φðT[figÞ φðTÞ � φðS[figÞ φðSÞ: This property establishes that the marginal gain player igets from joining a coalition is always greater when the coalition is larger. Therefore the Shapley values are always the greatest for the largest coalitions and that is why the method proposed is always doomed to mistake the largest subsets as communities. As we have pointed, the weakness is not on using cooperative games to define stable coalitions, but on using convex cooperative games. In this section, we will provide a simple and effective way to adjust this weakness. Our proposal is based on determining stability using a non-convex cooperative game. The computation of the expected weight on an arc. As we discussed in the previous section, weighted graph community games in which arc weights W ij �0 are convex games, so that they imply increasing values of the Shapley values and the tendency of detecting only large size communities. A straightforward way of avoiding convexity is considering an alternative set of weights, non necessarily non-negative, so that optimal stable coalitions of small size may emerge as well. Here, we propose to combine the weights defined by (4) with modularity, so that weights are normalized by their expected values and may take both negative and positive values. As a consequence, the resulting game is non-convex. The modularity function, see [1], is a well-known index to detect communities in networks. The index compares the edge density of the empirical graph G= (V,E) (unweighted and undirected), |E| = m, with the expected edge density of a theoretical graph G0= (V,E0) in which there are no communities by assumption. The expected edge density of G0is calculated using a null hypothesis, e.g. an assumption about the edge distribution, that is called the configuration model, [21]. If the graph does not contain communities, then for any given two nodes iand j with edge degrees k i and k j , the expected number of edges between iand jis approximated by kikj 2m. Let A ij = 1 if (i,j)2E,A ij = 0 otherwise (so that A= [A ij ] is the adjacency matrix of G). Moreover, let Pbe a partition of Vand let δ(i,j) be the Kronecker delta: δ(i,j) = 1 if i,j2V belong to the same community, δ(i,j) = 0 otherwise. Then the modularity function of a partition Pis: mðPÞ ¼ 1 2mX i;j2V Aij kikj 2m � �dði;jÞ:ð20Þ In the case under study, weights are defined through expression (4), in which the adjacency between nodes iand jis weighted by the common neighbors. However, modularity can be defined for weighted graphs as well. In the summation terms Aij kikj 2m � �, entries A ij are replaced by weights W ij ,k i replaced by weight sum W i =∑ j W ij , and mreplaced by W=∑ (i,j)2E W ij , as described in [22]. In this way, modularity is still a function that compares the actual indices of an empiric graph with the expected indices of a random graph. Using modularity, we can define modularity game (V,φ) as a weighted graph community game in which the PLOS ONE Overlapping communities detection through weighted graph community games PLOS ONE | https://doi.org/10.1371/journal.pone.0283857 April 4, 2023 8 / 35
characteristic function φis defined as in (1), but with the following weights: W0 ij ¼Wij WiWj 2W:ð21Þ In this case, W0 ij can take both positive and negative values, so that the game resulting from the characteristic function (1) is non-convex. We elaborate this model further, by noting that the modular term (21) should represent the difference between the empiric value W ij and its expected value under the assumption that the graph does not contain any communities. Unfortunately, the term WiWj 2Wis only an approximation of the true expectation and this can cause unexpected biases. For example, when weights W ij correspond to the adjacency matrix A ij 2{0, 1}, the term kikj 2mis an estimate of the probability of an arc between iand j, but, if the graph is unbalanced, the term can be greater than 1, which results in a non-sense estimation of this probability. In our application, expression (4) contains specific terms about the graph structure, such as the arcs and the common neighbours between two nodes, and potentially the bias between the true expectation and its approximation can be large. For this reason, we made a special effort in calculating the exact equation of the expected values of expression (4) under the assumption that there are no community in the graph. In [21], the random occurrence of a graph with no communities is calculated through the configuration model. The configuration model can be interpreted as the process of making a random graph with no communities through the following operations. Every arc e= (i,j) of the empirical graph G= (V,E) is cut into two parts, say l 1 and l 2 , with l 1 incident to iand l 2 incident to j, called stubs. Next, two different stubs are selected randomly and paired. We say that, if l 1 and l 2 are such stubs, then (l 1 ,l 2 ) is a match, e.g. an arc of the random graph G0= (V, E0). The way in which G0is built implies that the adjacency degree k i remains unvaried for all i, but eventual communities are broken by random pairings of stubs. Note that, from construction, we can interpret any occurrence of G0as a matching of 2mstubs. The process is exemplified in Fig 2. Here we show how to compute exactly the expected values of expression (4) using the configuration model. Expected weights depend on the the partition ratio P ij and the neighbourhood ratio CN ij of the random graphs obtained from the configuration model. By construction, the partition ratio P ij of the random graph is the same as the one of the empiric graph, but the neighbourhood ratio CN ij is different. To calculate CN ij , we introduce some notation. Recall that k i is the adjacency degree of node iand assume that the graph has medges. Let P adjacency (k i ,k j ,m) be the probability that node iand jare connected by an arc, let P common neighbour (k i ,k j ,k r ,m) be the probability that i and jare arc connected with r, so thar ris a common neighbor, and let P triangle (k i ,k j ,k r ,m) be the probability that iand jare arc connected and are also connected with r, so that the three Fig 2. Configuration model example. https://doi.org/10.1371/journal.pone.0283857.g002 PLOS ONE Overlapping communities detection through weighted graph community games PLOS ONE | https://doi.org/10.1371/journal.pone.0283857 April 4, 2023 9 / 35
arcs form a triangle. The notation emphasizes that probabilities depend on adjacency degrees k i ,k j ,k r and the total number of edges m. In the following proposition, we will derive closed form expressions for the above probabilities. Proposition 1.Let i, j, r be three nodes with adjacency degrees k i , k j , k r , respectively. Then, in the random graph configuration model Padjacencyðki;kj;mÞ ¼ X minfki;kjg t¼1ð 1Þtþ1 ki t � � kj t � �t! Qt p¼1ð2mþ12pÞ;ð22Þ Pcommon neighbourðki;kj;kr;mÞ ¼ X minfkj;kr1g t¼1ð 1Þtþ1Padjacencyðki;krt;mtÞ kr t � � kj t � �t! Qt p¼1ð2mþ12pÞ;ð23Þ Ptriangleðki;kj;kr;mÞ ¼ X minfki1;kj1g t¼1ð 1Þtþ1Pcommon neighbourðkit;kjt;kr;mtÞ ki t � � kj t � �t! Qt p¼1ð2mþ12pÞ:ð24Þ Proof. Applying the configuration model to G= (V,E), we obtain two stubs l 1 and l 2 , adjacent to iand j, respectively, for every arc e(i,j)2E. Then, we select two stubs at random and pair them until a random graph G0is obtained. Note that, from construction, we can interpret any occurrence of G0as a matching of 2mstubs. Given i,j2V, let S i = {l i(1) ,. . .,l i(ki )} be the set of stubs adjacent to iand S j = {l j(1) ,. . .,l j(kj )} be the set of stubs adjacent to j. Assuming a set of 2melements, there are ð2mÞ! 2mm!¼Qm p¼1ð2mþ 12pÞdifferent matching, see [23]. Therefore, if two stubs l 1 2S i and l 2 2S j are matched, there are Qm1 p¼1ð2m12pÞdifferent matching with the stubs remaining, because there are still 2m−2 stubs to pair. Due to this, the probability that two stubs l 1 and l 2 are joined, connecting nodes iand j, is: Qm1 p¼1ð2m12pÞ Qm p¼1ð2mþ12pÞ¼1 2m1:ð25Þ Next, we introduce random variables: Xl1l2¼ 1;if the stubs l1and l2are matched; 0;otherwise: 8 > > > < > > > : Obviously, the probability of Xl1l2¼1is PðXl1l2¼1Þ ¼ 1 2m1, as stated in (25). We can express the number of edges between two nodes iand jas the sum: X l12SiX l22Sj Xl1l2: The above expression represents the sum of the variables Xl1l2whose indices are one stub adjacent to iand another stub adjacent to j. Thus, the expected number of edges between iand PLOS ONE Overlapping communities detection through weighted graph community games PLOS ONE | https://doi.org/10.1371/journal.pone.0283857 April 4, 2023 10 / 35
jis: X l12SiX l22Sj E½Xl1l2� ¼ X l12SiX l22Sj PðXl1l2¼1Þ ¼ kikj 2m1: Note that in the modularity function (20), this value is approximated by kikj 2m. As we explain before, the expected number of edges is different to the probability of adjacency. The adjacency between two nodes iand jis the condition that there is at least one arc between iand jand it can be expressed as the union of the events fo:Xl1l2ðoÞ ¼ 1gwith l 1 2 S i and l 2 2S j , for the sake of simplicity, we refer to this set of events as fXl1l2¼1g. So, the adjacency probability of two nodes iand jis: P�[ l12Si l22Sj Xl1l2¼1 n o�:ð26Þ Let St ij be the set of all the different subsets of S i ×S j with size jSt ijj ¼ t. Applying the inclusion-exclusion law for the probability of union of events to expression (26), it follows that: P�[ l12Si l22Sj fXl1l2¼1g�¼X kikj t¼1ð 1Þtþ1X S2St ij P�\ ðl1;l2Þ2SfXl1l2¼1g�:ð27Þ By construction of the random graph G0, observe that the intersection of tdifferent sets fXl1l2¼1g, representing the match between stubs l 1 and l 2 , is empty if the same stub, l 1 or l 2 , is repeated more than once in different matches. Therefore, for each t, the non empty sets Tðl1;l2Þ2SfXl1l2¼1gthat appears in (27) are matching with tmatches. As a consequence, the summation on tis bounded to min{k i ,k j }, because the intersection of more than min{k i ,k j } different sets must repeat some stubs and so, its intersection is empty. Moreover, applying the same argument to calculate the probability of joining two stubs (25), the probability of joining tstubs from S i with other tstubs from S j is: Qmt p¼1ð2mþ12t2pÞ Qm p¼1ð2mþ12pÞ¼1 Qt p¼1ð2mþ12pÞ: Finally, to derive expected vales, we need to calculate the number of different subsets from S i ×S j with a size equal to tthat do not repeat any stubs. We have to consider tstubs from S i and tfrom S j , and then all the possible matchings between stubs of different sets. There are ki t � different subsets of tstubs from S i and kj t ��different subsets of tstubs from S j . We can match the tstubs of one set with the other tstubs of the other set in t! different ways, obtaining the following expression for the probability of events ensuring that node iand jare connected, in PLOS ONE Overlapping communities detection through weighted graph community games PLOS ONE | https://doi.org/10.1371/journal.pone.0283857 April 4, 2023 11 / 35
overlapping communities reveals nodes that are structurally different from others, forming the bulk of a core/periphery separation. First, we apply the model F∗ ShMod to the the Highland tribes network, see [26]. First, model F∗ ShMod is run with p= 1 and results are in Fig 7a. There, it can be seen that, if no overlapping communities are allowed, then the model detects one community composed of all the nodes. Conversely, model F∗ ShMod is run with parameters n c = 3 and p= 2, results are reported in Fig 7b. It can be seen that the role of different nodes is emerged. There, three communities of different size have been detected, with some nodes (the red ones) belonging to more than one community forming the core of the system of alliances. Next, we apply model F∗ ShMod to the Windsurfers network, see [27]. Run with parameter p= 1, the model detected the two communities reported in Fig 8a. Run with parameters n c = 2 and p= 2, the model detected the communities reported in Fig 8b. As can be seen, the results with overlapping communities are a refinement of the disjoint communities. Nodes that are in the border between the two groups are highlighted as members of both, forming the bulk of a core/periphery network segmentation. To summarize our findings, the test of models F∗ ShMod on four typical benchmark networks revealed: • Results between F∗ ShMod and F0 ShMod are different. As the latter is an approximation of the former, it reveals that the contribution of Theorem 1 to model development is substantial. • Results between non-overlapping and overlapping community models are different. The former can reveal not only group membership, but nodes that could act as potential bridges between communities. Fig 7. Highland tribes community structures. Community structures obtained by (a) F∗ ShMod with parameters n c =n, p= 1, (b) F∗ ShMod with parameters n c = 3, p= 2. https://doi.org/10.1371/journal.pone.0283857.g007 Fig 8. Windsurfers community structures. Community structures obtained by (a) F∗ ShMod with parameters n c =n, p= 1, (b) F∗ ShMod with parameters n c = 2, p= 2. https://doi.org/10.1371/journal.pone.0283857.g008 PLOS ONE Overlapping communities detection through weighted graph community games PLOS ONE | https://doi.org/10.1371/journal.pone.0283857 April 4, 2023 18 / 35
Local Stability Exploration: An heuristic algorithms to detect overlapping communities Problems F∗ ShMod and F0 ShMod are Integer Linear Programming (ILP) models whose solution computational times can be impractical when the instances to solve are large. This is normal when we deal with a NP-hard problem as the case of communities detection. Nevertheless, for large instances the ILP formulation can be applied to devise heuristic algorithms that could approximate the optimal solution in short computing time. Here we propose a method, that we will call Local Stability Exploration (LSE), that is based on local search. Suppose that a set of feasible communities P¼ fS1;...;Sncg;Si�V;i¼1;...;ncis given, we will call such Pan incumbent solution. Pfeasible means that it satisfies the ILP model constraints, so that i) every node belongs to at least one community, Snc k¼1Sk¼V, ii) there is not strict inclusion between communities, ∄k;r¼1;. . . ;nc;k6¼ r;such that S k �S r , iii) the maximum number of communities to which a node can belong is not exceeded by any node, i.e. 8i2Vthe inequality |{k= 1, . . .,n c :i2S k }| �pis fulfilled; and iv) all communities are stable. Next, we try to modify Pto obtain a new feasible solution P0with an improved objective function. We consider three possible modification of P, obtained by moves that are called Add, Remove, and Swap. Add is the move that joins a node to a community, allowing in this way multiple communities assignments. Remove is the move that takes away a node from a community. Swap is the move that switch two nodes between two communities. These moves are applied if and only if the new obtained P0is feasible. That is, after a move it must not occur that 1) a node does not belong to any community 2) a node belongs to more communities than allowed, maximum number of communities pto which a node can belong; 3) one community is included in another, 4) modified communities are not stable. For a feasible starting solution, the procedure is summarized in Algorithm 1. There, the triplet (i,k, 1) is the move of adding node ito community k, the triplet (i,k, 2) is the move of removing node ifrom community k, the 5-tuple (i,k,i0,k0, 3) is swapping nodes iand i0 between communities kand k0. It can be seen that from Line 9 to Line 22 all feasible moves are considered. In Lines 12, 15 and 20 the increases of the objective function are calculated using the following notation: Let C i = {k2{1, . . .,n c }:i2S k }, that is, C i is the index set of the communities to which ibelongs, then the objective function can be written as: f∗Pð Þ ¼ X n i;j2V i<j Ci\Cj6¼ ; W∗ ij Note that the condition C i \C j 6¼;is the condition that there is at least one community to which both iand jbelong to. However, from the computational efficiency it is better to calculate just the increase of the objective function, as is done in lines 12, 15, 20. The new solution P0is the one that obtains the maximum increase. The algorithm stops when condition of Line 42 applies, as there are no improvements and a local optimum has been reached. Algorithm 1 Local stability exploration algorithm 1: procedure LOCAL STABILITY EXPLORATION 2: P¼ fS1;. . . ;Sncg Initial Stable Communities ⊳Π is obtained by peculiar subroutines 3: for iin Vdo 4: C i = {k2{1, . . .,n c }: i2S k } 5: end for 6: f P n i;j2V i<j Ci\Cj6¼ ; W∗ ij ⊳Objective function PLOS ONE Overlapping communities detection through weighted graph community games PLOS ONE | https://doi.org/10.1371/journal.pone.0283857 April 4, 2023 19 / 35
7: local_opt =FALSE ⊳Condition for a local optimum 8: while local_opt =FALSE do 9: Δ Feasible_Moves(Π) ⊳Δ: list of admissible moves for Π. 10: for (i,k,d) in Δ do 11: if d = 1 then 12: dikd P j2Sk Ci\Cj¼ ; W∗ ij 13: end if 14: if d = 2 then 15: dikd P j2Sknfig jCi\Cjj ¼ 1 W∗ ij 16: end if 18: end for 18: for (i,k,i0,k0, 3) 2Δdo 19: if d=3then 20: diki0k0d P j2Sk0nfi0g Ci\Cj¼ ; W∗ ij P j2SknSk0Sif gð Þ jCi\Cjj ¼ 1 W∗ ij þP j2Sknif g Ci0\Cj¼ ; W∗ i0jP j2Sk0nSkSi0 f gð Þ jCi0\Cjj ¼ 1 W∗ i0j 21: end if 22: end for 23: (i*,k*,d*)2argmax{δ ikd |(i,k,d)2Δ} ⊳Select the move that increases the most 24: (i*,k*,i0*,k0*,d*)2argmax{δ iki 0 k 0 d |(i,k,i0,k0,d)2Δ} ⊳Select the move that increases the most 25: if δ i*k*i 0 *k 0 *d* > max{0, δ i*k*d* }then 26: f f+δ i*k*i 0 *k 0 *d* ⊳Update f 27: S k* S k* [{i0*}\{i*} 28: S k 0 * S k 0 * [{i*}\{i0*}⊳Update Π 29: C i* C i* [{k0*}\{k*} 30: C i 0 * C i 0 * [{k*}\{k0*} 31: else 32: if δ i*k*d* > 0 then 33: f f+δ i*k*d* ⊳Update f 34: if d*= 1 then 35: S k* S k* [{i*}⊳Update Π 36: C i* C i* [{k*} 37: else 38: S k* S k* \{i*}⊳Update 39: C i* C i* \{k*} 40: end if 41: else 42: local_opt =TRUE 43: end if 44: end if 45: end while 46: return Π⊳Return the local optimum 47: end procedure It remains to comment how feasible starting solutions can be obtained in Line 2 of Algorithm LSE. Depending on problems, we tested various procedures. The first possibility is to start with an unfeasible solution P, because it contains unstable communities. Then Algorithm LSE is run without imposing that new solutions P0should be stable, but once that a feasible one has been found, then all forthcoming solutions must remain feasible too. The first unfeasible Pcan be a random assignment to communities, but another possibility is solving F∗ ShMod for p= 1, that is, when overlapping is not allowed, as the problem is usually solved faster than the cases in which p>1. Another possibility that has been used for the problems with the PLOS ONE Overlapping communities detection through weighted graph community games PLOS ONE | https://doi.org/10.1371/journal.pone.0283857 April 4, 2023 20 / 35
largest size is solving F∗ ShMod by branch-and-bound, but stop the search when the first feasible solution has been found and next using it as the starting solution in Line 2. All methods can be combined using any multi-start strategy, that is, repeating Algorithm 1 many times with different starting solutions to obtain sufficient diversification and exploration of the solution space. Finally, Algorithm LSE has been explained to solve model F∗ ShMod, but it can be applied to F0 ShMod with straightforward modifications. A preliminary test of the quality of the LSE algorithm has been run on the previous networks. We run a multi-start version allowing t max = 10 starting solutions each run. Results about computational times and solution quality for different parameters configurations are reported in Table 1. It can be seen that the LSE heuristic algorithm reduces the computing time significantly with respect to the ILP solution for both models F∗ ShMod and F0 ShMod, while the optimal solution has been achieved in all the cases but one. Moreover, we applied the LSE algorithm to some large-scale real data sets that are impractical for any ILP model, in order to test the scalability of our heuristic. The solved data sets are the American college football network with 115 nodes, see [28], the Jazz musician network with 198 nodes, see [29], and C. metabolic network with 453 nodes, see [30]. These real data sets examples are commonly used in literature. The exact expected weights W∗ ij cannot be computed for graphs with a large number of edges, so we used the approximated expected weights W0 ij. We report the results of our methods in Figs 9–11. Results and discussion We are going to analyze the main features of the ILP models F∗ ShMod,F0 ShMod and the heuristic Algorithm 1 when they are applied to medium and large size networks, most precisely, Table 1. Computational results of the solution methods. Dataset n c pModel Solving method Time (s) Objective value Zachary’s karate club 4 2 F∗ ShMod Exact ILP 1530 162.469 LSE heuristic 77 162.469 Zachary’s karate club 4 2 F0 ShMod Exact ILP 316 129.39 LSE heuristic 6 129.279 Zachary’s karate club 3 2 F∗ ShMod Exact ILP 139 157.652 LSE heuristic 15 157.652 Zachary’s karate club 3 2 F0 ShMod Exact ILP 53 122.578 LSE heuristic 6 122.578 Highland tribes 3 2 F∗ ShMod Exact ILP 31 89.654 LSE heuristic 4 89.654 Highland tribes 3 2 F0 ShMod Exact ILP 8 47.4516 LSE heuristic 0.46 47.4516 Zebra communication 3 2 F∗ ShMod Exact ILP 12 313.869 LSE heuristic 2.8 313.869 Zebra communication 3 2 F0 ShMod Exact ILP 16 146.858 LSE heuristic 2.2 146.858 Windsurfers 2 2 F∗ ShMod Exact ILP 1059 583.388 LSE heuristic 78 583.388 Windsurfers 2 2 F0 ShMod Exact ILP 14 292.662 LSE heuristic 7.5 292.662 https://doi.org/10.1371/journal.pone.0283857.t001 PLOS ONE Overlapping communities detection through weighted graph community games PLOS ONE | https://doi.org/10.1371/journal.pone.0283857 April 4, 2023 21 / 35
whether they can detect the true overlapping communities of randomly generated networks, as it is done in [31]. Random networks are generated using the procedure proposed in [16], but with some variations to allow for communities that overlap. Most peculiarly, in our simulation we must distinguish between bridge and non-bridge nodes, the former being the nodes that belongs to more than one community. The main parameters characterizing the simulated networks are: •N: the number of nodes. •n c : the number of communities. Fig 10. Jazz music community structure. Community structure obtained by LSE heuristic with W0 ij weights and parameters n c = 6, p= 2. https://doi.org/10.1371/journal.pone.0283857.g010 Fig 9. American college football community structure. Community structure obtained by LSE heuristic with W0 ij weights and parameters n c = 7, p= 2. https://doi.org/10.1371/journal.pone.0283857.g009 PLOS ONE Overlapping communities detection through weighted graph community games PLOS ONE | https://doi.org/10.1371/journal.pone.0283857 April 4, 2023 22 / 35
•p: the maximum number of communities to which a node can belong to. •N o : the number of nodes that belongs to more than one communities, that is, they are bridges. Next, communities are defined by the probability by which community nodes can establish a link between themselves. Those probabilities are controlled by parameters: • 1 −μ: fraction of links between non-bridge nodes belonging to the same community. • 1 −μ o : fraction of links between bridge nodes and other nodes of the communities where the bridge node belongs to. There are other parameters characterizing the simulated networks, such as the number of arcs, the node degrees, the community sizes and so on, whose purpose is to simulate networks with the same characteristics of the empiric ones. We report all these features in S1 Appendix, with the pseudo-code describing our implementation of Lancichenetti et al. algorithm. The solution quality of our models is measured comparing their results with the true community structures (known by simulation). True and estimated structure may differ for: • The community composition; • The identification of the bridge nodes. The statistics to compare the community composition are: • the Normalized Mutual Information (NMI) index for overlapping partitions, presented in [32]; • the Omega index (OI), presented in [33]. Both statistics range between 0 and 1, with values closer to 1 indicating strong correspondence between true and estimated communities. Fig 11. C. metabolic community structure. Community structure obtained by LSE heuristic with W0 ij weights and parameters n c = 10, p= 2. https://doi.org/10.1371/journal.pone.0283857.g011 PLOS ONE Overlapping communities detection through weighted graph community games PLOS ONE | https://doi.org/10.1371/journal.pone.0283857 April 4, 2023 23 / 35
The statistics to compare the identification of bridge nodes are based on a set of indices which depend on the values of the confusion matrix associated to the identification of bridge nodes. Each element of the confusion matrix is defined as follows • True Positive (TP): Nodes successfully detected as bridge. • True Negative (TN): Nodes successfully detected as non-bridge. • False Positive (FP): Nodes wrongly detected as bridge. • False Negative (FN): Nodes wrongly detected as non-bridge. Then, we consider the following indices. the accuracy defined as TPþTN TPþTNþFPþFN, the True Positive Rate (TPR):TPR ¼TP TPþFN, the False Positive Rate (FPR):FPR ¼FP TNþFP, the Area Under Curve (AUC):AUC ¼1FPRþTPR 2, the Precision defined as TP TPþFP, the F1 score:F1¼2TP 2TPþFPþFN; Test 1: Detecting non overlapping communities: As a first test, we apply the ILP models F∗ ShMod,F0 ShMod and the Algorithms LSE to the case in which communities do not overlap, that is, p= 1, to see whether the approximate result of algorithm LSE are reliable, with respect to what is found by the respective optimal ILP models. The ILP solution of F∗ ShMod and F0 ShMod can be obtained in short computational times only for moderate size networks, so we consider N= 40, 60 to solve within the time limit of 100 or 200 seconds respectively. The LSE heuristic has been run with t max = 5 multiple starting solution, guaranteeing that its computational times are a fraction of the exact method. For fixed Nand n c , we let μ= 0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, as in [16] to control for the effect of mixing parameter. For each parameter set, either 50 or 100 random networks are generated and indices are calculated as averages on all instances. Results are reported in Table 2. The first two rows of this table give the ILP formulation (F∗ ShMod or F0 ShMod) used in the corresponding method: exact (ILP) or (LSE) heuristic to provide an initial solution. The third row describes the parameters of the instances (N,n c ,μ) and the index reported below (NMI or Omega). By columns, the layout of this table is organized in three blocks. The first one with three columns describes the instances. The next two blocks, each one with four columns, report the average values of the NMI and Omega indices for each combination of solution method. Results in bold report the best behaviour among similar index for the corresponding solution methods. One can easily observe that using formulation F∗ ShMod in the ILP or in the LSE heuristic provides better solutions than F0 ShMod. For each combinations of parameters Nand n c , the NMI and OI of each solution method are also shown as a function of μin Figs 12–14 to compare the formulations F∗ ShMod and F0 ShMod. The exact formulation F∗ ShMod obtains, in general, better NMI results and also better OI results in more cases than F0 ShMod; except for N= 40 and n c = 4. In this case, the behaviour of OI is similar in both formulations. However, also for N= 40, the exact solution of model F∗ ShMod is superior to the other two approaches, namely the heuristic LSE and the exact model F0 ShMod, as the curves of the NMI and Omega statistics are above the others for most values of N,n c and μ. PLOS ONE Overlapping communities detection through weighted graph community games PLOS ONE | https://doi.org/10.1371/journal.pone.0283857 April 4, 2023 24 / 35
When μis above the threshold 0.3, the solution quality of the method deteriorates for the joint effect of two factors: 1) communities are less well-separated, 2) exact solution has not been obtained within the considered time limit. However, this is not an actual drawback since for those parameter values, communities are essentially meaningless. Test 2: Detecting overlapping communities on small networks: Networks with overlapping communities have been simulated with the same parameters used before, but now communities can overlap. We control the overlap with parameters p= {2, 3} and μ o = {0.5, 0.7}. The choice of these parameters is justified since for p= 2 the smallest possible μ o value is 0.5 and for p= 3 the smallest possible μ o value is approximately 0.7. Moreover, the number of bridge Table 2. Computational results about networks with non-overlapping communities. Model F∗ ShMod F0 ShMod Method ILP LSE ILP LSE Nn c μNMI OI NMI OI NMI OI NMI OI 40 6 0 0.95 0.99 0.88 0.91 0.95 1 0.82 0.85 0.1 0.96 0.98 0.88 0.94 0.95 10.82 0.86 0.2 0.91 0.84 0.88 0.91 0.79 0.86 0.79 0.83 0.3 0.85 0.79 0.82 0.85 0.66 0.72 0.73 0.78 0.4 0.62 0.47 0.7 0.7 0.45 0.5 0.59 0.63 0.5 0.38 0.32 0.48 0.5 0.16 0.19 0.46 0.49 0.6 0.39 0.2 0.29 0.3 0.1 0.11 0.28 0.3 40 4 0 0.88 0.99 0.87 0.98 0.88 1 0.83 0.92 0.1 0.88 0.99 0.86 0.99 0.88 1 0.81 0.91 0.2 0.88 0.94 0.85 0.97 0.88 0.99 0.83 0.92 0.3 0.82 0.79 0.83 0.93 0.78 0.88 0.81 0.9 0.4 0.63 0.63 0.78 0.87 0.57 0.68 0.76 0.85 0.5 0.35 0.34 0.6 0.71 0.17 0.21 0.57 0.65 0.6 0.15 0.23 0.36 0.41 0.07 0.09 0.34 0.41 60 6 0 0.79 0.85 0.89 0.98 0.79 0.87 0.83 0.92 0.1 0.72 0.78 0.9 0.99 0.53 0.6 0.84 0.93 0.2 0.7 0.73 0.9 0.98 0.33 0.41 0.83 0.92 0.3 0.63 0.65 0.89 0.97 0.07 0.08 0.82 0.9 0.4 0.36 0.39 0.87 0.93 0 0.01 0.79 0.87 0.5 0.28 0.3 0.73 0.79 0 0 0.66 0.74 0.6 0.13 0.15 0.44 0.52 0 0 0.39 0.47 Maximum NMI and OI values for each combination of parameters and each method are highlighted in bold. https://doi.org/10.1371/journal.pone.0283857.t002 Fig 12. Test results on non-overlapping communities, parameters N= 40, n c = 6. (a) Average NMI for each solution method, (b) Average Omega for each solution method. https://doi.org/10.1371/journal.pone.0283857.g012 PLOS ONE Overlapping communities detection through weighted graph community games PLOS ONE | https://doi.org/10.1371/journal.pone.0283857 April 4, 2023 25 / 35
nodes N o is approximately 10% of all the nodes, and we change this value to asses how it affects the computational results. Problems with overlapping communities are harder to solve, therefore we limit the graph size to N= 40 and increase the time limit to 200 seconds. Table 3 reports the computational results with a layout similar to Table 2. It can be seen that the best values of both the Omega and NME indices are obtained with the LSE heuristic, applied to the F0 ShMod formulation. The LSE heuristic applied to F∗ ShMod provides the second best results (with a few exceptions in which it becomes the best one) and the third one is the ILP formulation. The reason of the poor performance of the ILP methods is due to the fact that they were not able to terminate the computation in the imposed time limit and the solution that they provide is far from optimality. Results of Table 3 are graphically reported in Figs 15–19, where it can be seen that the purple and green curve, representing the LSE heuristics, are very close with each other and they are much above the result of the truncated ILP. It is also noteworthy that weights from F∗ ShMod improve the results of the ILP method. Fig 14. Test results on non-overlapping communities, parameters N= 60, n c = 6. (a) Average NMI for each solution method, (b) Average Omega for each solution method. https://doi.org/10.1371/journal.pone.0283857.g014 Fig 13. Test results on non-overlapping communities, parameters N= 40, n c = 4. (a) Average NMI for each solution method, (b) Average Omega for each solution method. https://doi.org/10.1371/journal.pone.0283857.g013 PLOS ONE Overlapping communities detection through weighted graph community games PLOS ONE | https://doi.org/10.1371/journal.pone.0283857 April 4, 2023 26 / 35
Test 3: Detecting overlapping communities on large-scale networks: In the last experiment, we have applied the LSE heuristics, using both the F∗ ShMod and F0 ShMod models, to the largest networks composed of 500 or 1000 nodes. As before, we control the overlap between communities with parameters p= {2, 3} and μ o = {0.6, 0.7}, the number of bridge nodes are N o = {20, 50}. In Table 4, we report the NMI and OI statistics calculated by the two methods. It can be seen that they have lower values than what obtained in the smallest networks, due the fact that Table 3. Computational results about networks with overlapping communities. Model F∗ ShMod F0 ShMod Method ILP LSE ILP LSE n c p μ o N o μNMI OI NMI OI NMI OI NMI OI 4 2 0.5 1 0 0.88 0.92 0.86 0.92 0.88 0.87 0.91 0.95 0.1 0.76 0.78 0.88 0.91 0.71 0.73 0.9 0.95 0.2 0.63 0.66 0.85 0.91 0.45 0.46 0.88 0.92 0.3 0.57 0.63 0.82 0.88 0.28 0.3 0.83 0.88 0.4 0.4 0.46 0.76 0.81 0.17 0.18 0.72 0.78 0.5 0.28 0.33 0.49 0.57 0.1 0.12 0.5 0.57 0.6 0.13 0.16 0.29 0.35 0.05 0.05 0.31 0.37 3 0 0.88 0.88 0.9 0.92 0.81 0.82 0.95 0.96 0.1 0.77 0.78 0.91 0.93 0.66 0.67 0.95 0.96 0.2 0.61 0.72 0.88 0.91 0.5 0.53 0.9 0.93 0.3 0.55 0.57 0.85 0.89 0.32 0.33 0.82 0.85 0.4 0.52 0.46 0.78 0.82 0.19 0.21 0.72 0.75 0.5 0.26 0.29 0.46 0.53 0.06 0.08 0.5 0.55 0.6 0.17 0.19 0.32 0.38 0.07 0.08 0.33 0.38 5 0 0.74 0.79 0.91 0.93 0.72 0.71 0.94 0.95 0.1 0.67 0.69 0.9 0.91 0.66 0.67 0.96 0.97 0.2 0.58 0.65 0.91 0.92 0.42 0.46 0.9 0.91 0.3 0.53 0.56 0.85 0.87 0.28 0.31 0.82 0.85 0.4 0.42 0.42 0.7 0.76 0.15 0.17 0.64 0.68 0.5 0.28 0.27 0.44 0.52 0.06 0.06 0.44 0.49 0.6 0.15 0.21 0.31 0.37 0.06 0.07 0.28 0.32 0.7 3 0 0.87 0.8 0.88 0.9 0.94 0.94 0.94 0.96 0.1 0.77 0.68 0.88 0.9 0.79 0.8 0.95 0.96 0.2 0.71 0.7 0.88 0.91 0.46 0.48 0.89 0.91 0.3 0.55 0.54 0.83 0.87 0.24 0.26 0.85 0.88 0.4 0.42 0.38 0.72 0.77 0.18 0.21 0.67 0.71 0.5 0.23 0.29 0.49 0.56 0.08 0.09 0.47 0.51 0.6 0.17 0.15 0.33 0.39 0.08 0.09 0.29 0.34 3 0.7 3 0 0.74 0.7 0.82 0.84 0.8 0.83 0.88 0.89 0.1 0.67 0.64 0.85 0.87 0.57 0.6 0.87 0.89 0.2 0.59 0.56 0.83 0.86 0.3 0.34 0.83 0.85 0.3 0.5 0.42 0.76 0.8 0.18 0.2 0.71 0.76 0.4 0.32 0.41 0.63 0.69 0.17 0.21 0.56 0.62 0.5 0.22 0.26 0.35 0.43 0.11 0.12 0.37 0.44 0.6 0.14 0.16 0.21 0.27 0.07 0.08 0.22 0.26 Maximum NMI and OI values for each combination of parameters and each method are highlighted in bold. https://doi.org/10.1371/journal.pone.0283857.t003 PLOS ONE Overlapping communities detection through weighted graph community games PLOS ONE | https://doi.org/10.1371/journal.pone.0283857 April 4, 2023 27 / 35
Methodology: Stefano Benati, Justo Puerto, Antonio M. Rodrı ´guez-Chı ´a, Francisco Temprano. Project administration: Stefano Benati, Justo Puerto, Antonio M. Rodrı ´guez-Chı ´a, Francisco Temprano. Resources: Stefano Benati, Justo Puerto, Antonio M. Rodrı ´guez-Chı ´a, Francisco Temprano. Software: Stefano Benati, Justo Puerto, Antonio M. Rodrı ´guez-Chı ´a, Francisco Temprano. Supervision: Stefano Benati, Justo Puerto, Antonio M. Rodrı ´guez-Chı ´a. Validation: Stefano Benati, Justo Puerto, Antonio M. Rodrı ´guez-Chı ´a, Francisco Temprano. Visualization: Stefano Benati, Justo Puerto, Antonio M. Rodrı ´guez-Chı ´a, Francisco Temprano. Writing – original draft: Stefano Benati, Justo Puerto, Antonio M. Rodrı ´guez-Chı ´a, Francisco Temprano. Writing – review & editing: Stefano Benati, Justo Puerto, Antonio M. Rodrı ´guez-Chı ´a, Francisco Temprano. References 1. Girvan M, Newman MEJ. Finding and evaluating community structure in networks. Phys Rev E. 2004;( 69 (2), 026113). https://doi.org/10.1103/PhysRevE.69.026113 PMID: 14995526 2. Fortunato S, Hric D. Community detection in networks: A user guide. Physics Reports. 2016; 659:1–44. https://doi.org/10.1016/j.physrep.2016.09.002 3. Palla G, Dere ´nyi I, Farkas I, Vicsek T. Uncovering the Overlapping Community Structure of Complex Networks in Nature and Society. Nature. 2005;( 435 (7043)). https://doi.org/10.1038/nature03607 PMID: 15944704 4. Xie J, Kelley S, Szymanski BK. Overlapping Community Detection in Networks: The State-of-the-Art and Comparative Study. Comput Surv. 2013; 45(43):1–35. https://doi.org/10.1145/2501654.2501657 5. Agarwal G, Kempe D. Modularity-maximizing graph communities via mathematical programming. The European Physical Journal B. 2008; 66(3):409–418. https://doi.org/10.1140/epjb/e2008-00425-1 6. Li Z, Zhang XS, Wang RS, Liu H, Zhang S. Discovering Link Communities in Complex Networks by an Integer Programming Model and a Genetic Algorithm. PLoS ONE. 2013;( 8 (12), e83739). https://doi. org/10.1371/journal.pone.0083739 PMID: 24386268 7. Bennett L, Kittas A, Liu S, Papageorgiou LG, Tsoka S. Community Structure Detection for Overlapping Modules through Mathematical Programming in Protein Interaction Networks. PLoS ONE. 2014;( 9(11): e112821). https://doi.org/10.1371/journal.pone.0112821 PMID: 25412367 8. Costa A, Ng TS, Foo LX. Complete mixed integer linear programming formulations for modularity density based clustering. Discrete Optimization. 2017;(25):141–158. https://doi.org/10.1016/j.disopt.2017. 03.002 9. Zhang S, Wang RS, Zhang X. Identification of overlapping community structure in complex networks using fuzzy c-means clustering. Physica A. 2007;(374):483–490. https://doi.org/10.1016/j.physa.2006. 07.023 10. Nepusz T, Petroczi A, Negyessy L, Bazso F. Fuzzy Communities and the Concept of Bridgeness in Complex Networks. Physical Review E. 2008; 77:16–107. https://doi.org/10.1103/PhysRevE.77. 016107 PMID: 18351915 11. Nicosia V, Mangioni G, Carchiolo V, Malgeri M. Extending the definition of modularity to directed graphs with overlapping communities. J Stat Mech Theory Exp. 2009;((03) (2009) P03024). 12. Chen D, Shang M, Fu Y. Detecting overlapping communities of weighted networks via a local algorithm. Physica A: Statistical Mechanics and its Applications. 2010;(389):4177–4187. https://doi.org/10.1016/j. physa.2010.05.046 13. Chitra Devi J, Poovammal E. An Analysis of Overlapping Community Detection Algorithms in Social Networks. Procedia Computer Science. 2016;(89):349–358. https://doi.org/10.1016/j.procs.2016.06. 082 PLOS ONE Overlapping communities detection through weighted graph community games PLOS ONE | https://doi.org/10.1371/journal.pone.0283857 April 4, 2023 34 / 35
14. Benati S, Puerto J, Rodrı ´guez-Chı ´a AM, Temprano F. A mathematical programming approach to overlapping community detection. Physica A: Statistical Mechanics and its Applications. 2022; 602:127628. https://doi.org/10.1016/j.physa.2022.127628 15. Jonnalagadda A, Kuppusamy L. A cooperative game framework for detecting overlapping communities in social networks. Physica A. 2018;(491):498–515. https://doi.org/10.1016/j.physa.2017.08.111 16. Lancichinetti A, Fortunato S, Radicchi F. Benchmark graphs for testing community detection algorithms. Phys Rev E. 2008; 78:046110. https://doi.org/10.1103/PhysRevE.78.046110 PMID: 18999496 17. Demange G. Intermediate preferences and stable coalition structures. Journal of Mathematical Economics. 1994; 23(1):45–58. https://doi.org/10.1016/0304-4068(94)90035-3 18. Carraro C, Marchiori C. DP3258 Stable Coalitions. CEPR Press Discussion Paper. 2002; (3258). 19. D’Aspremont C, Jacquemin A, Gabszewicz JJ, Weymark JA. On the Stability of Collusive Price Leadership. The Canadian Journal of Economics / Revue canadienne d’Economique. 1983; 16(1):17–25. https://doi.org/10.2307/134972 20. Caparros A, Giraud-He ´raud E, Hammoudi A, Tazdaït T. Coalition Stability with Heterogeneous Agents. Economics Bulletin. 2011; 31(1):286–296. 21. Newman MEJ. Networks: an introduction. Oxford University Press; 2010. 22. Newman MEJ. Analysis of weighted networks. Phys Rev E. 2004;( 70 (5), 056131). https://doi.org/10. 1103/PhysRevE.70.056131 PMID: 15600716 23. Callan D. A combinatorial survey of identities for the double factorial. 2009;. 24. Zachary WW. An Information Flow Model for Conflict and Fission in Small Groups. Journal of Anthropological Research. 1977; 33(4):452–473. https://doi.org/10.1086/jar.33.4.3629752 25. Sundaresan SR, Fischhoff IR, Dushoff J, Rubenstein DI. Network metrics reveal diVerences in social organization between two Wssion-fusion species, Grevy’s zebra and onager. Oecologia. 2007; 151:140–149. https://doi.org/10.1007/s00442-006-0553-6 PMID: 16964497 26. Read KE. Cultures of the central highlands, New Guinea. Southwestern Journal of Anthropology. 1954; p. 1–43. https://doi.org/10.1086/soutjanth.10.1.3629074 27. Freeman LC, Freeman SC, Michaelson AG. On human social intelligence. Journal of Social Biological Structure. 1988; 11:415–425. https://doi.org/10.1016/0140-1750(88)90080-2 28. Girvan M, Newman MEJ. Community structure in social and biological networks. Proceedings of the National Academy of Sciences. 2002; 99(12):7821–7826. https://doi.org/10.1073/pnas.122653799 PMID: 12060727 29. Gleiser P, Danon L. Community Structure in Jazz. Advances in Complex Systems (ACS). 2003; 06:565–573. https://doi.org/10.1142/S0219525903001067 30. Jeong H, Tombor B, Albert R, Oltvai Z, Barabasi AL. The Large-Scale Organization of Metabolic Networks. Nature. 2000; 407(6804):651–654. https://doi.org/10.1038/35036627 PMID: 11034217 31. Tandon A, Albeshri A, Thayananthan V, Alhalabi W, Radicchi F, Fortunato S. Community detection in networks using graph embeddings. Phys Rev E. 2021; 103:022316. https://doi.org/10.1103/PhysRevE. 103.022316 PMID: 33736102 32. Lancichinetti A, Fortunato S, Kerte ´sz J. Detecting the overlapping and hierarchical community structure in complex networks. New Journal of Physics. 2009; 11(3). https://doi.org/10.1088/1367-2630/11/3/ 033015 33. Collins LM, Dent CW. Omega: A General Formulation of the Rand Index of Cluster Recovery Suitable for Non-disjoint Solutions. Multivariate Behavioral Research. 2005; 23:231–242. https://doi.org/10. 1207/s15327906mbr2302_6 PLOS ONE Overlapping communities detection through weighted graph community games PLOS ONE | https://doi.org/10.1371/journal.pone.0283857 April 4, 2023 35 / 35
Chapter 8 Mixed-integer linear programming formulations and column generation algorithms for the Minimum Normalized Cuts problem on networks 99
European Journal of Operational Research 316 (2024) 519–538 Available online 1 March 2024 0377-2217/© 2024 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). Contents lists available at ScienceDirect European Journal of Operational Research journal homepage: www.elsevier.com/locate/eor Discrete Optimization Mixed-integer linear programming formulations and column generation algorithms for the Minimum Normalized Cuts problem on networks Diego Ponce, Justo Puerto, Francisco Temprano∗ IMUS and Stats & OR Department, Universidad de Sevilla, Spain ARTICLE INFO MSC: 05C82 90C27 90C30 90C35 91C20 Keywords: Normalized cuts Complex networks Community detection Mathematical programming Image segmentation ABSTRACT This paper deals with the 𝑘-way normalized cut problem in complex networks. It presents a methodology that uses mathematical optimization to provide mixed-integer linear programming formulations for the problem. The paper also develops a branch-and-price algorithm for the above-mentioned problem which scales better than the compact formulations. Additionally, a heuristic algorithm which is able to approximate largescale image problems in those cases where the exact methods are not applicable is presented. Extensive computational experiments assess the usefulness of these methods to solve the 𝑘-way normalized cut problem. Finally, we have applied the minimum normalized cut objective function to the segmentation of actual images, showing the applicability of the introduced methodology. 1. Introduction The problem of partitioning a graph into a given number of 𝑘connected components is a long studied topic in the literature (Hansen & Jaumard,1997). It arises in a wide variety of practical applications in different contexts, such as engineering, economics, logistics, psychology, medicine, and electoral systems (Arredondo et al.,2021;Maravalle et al.,1997;Matić & Grbić,2020;Ricca et al.,2013;Validi et al., 2022;Zhou et al.,2019). One of its most interesting applications is image segmentation which can be obtained by grouping elements of images by key factors such as similarity, proximity, and good continuation (Wertheimer,1938). Nowadays, there are different approaches to the image segmentation problem based on several clustering techniques (among others), but there not yet a common agreement on what is a correct or optimal solution. One of the most used and well-recognized techniques in this field is based on the minimum cut problem. The Minimum Cut Problem, in graph theory, provides a bipartition that minimizes the number of edges between nodes from different subsets. This classical problem can be extended to the minimum 𝑘-cut problem (Goldschmidt & Hochbaum, 1994) that minimizes the number of edges between nodes from different subsets of a 𝑘-partition. Furthermore, it can be extended to weighted and directed graphs. ∗Corresponding author. E-mail addresses: [email protected] (D. Ponce), [email protected] (J. Puerto), [email protected] (F. Temprano). The advantage of the minimum cut problem is that it can be used as a generic clustering method when the data to be partitioned are codified in a network or another type of structure that can be transformed into a graph. Data samples from a metric space can be represented as nodes and weighted edges whose weights are inversely proportional to the distances between them. Thus, minimum cut problem is considered a community detection method that attempts to partition a network into disjoint groups densely connected among themselves and poorly connected to other communities. Nevertheless, the main drawback of a minimum cut approach is that the internal number of edges of the communities is not taken into account, so this approach tends to detect small communities or communities that contain nodes with a small adjacency degree. This is due to the reduced number of edges between communities with these properties. Therefore, the minimum cut sometimes leads to wrong community structures, such as in the following example. Example 1. The bipartition that represents the minimum cut of Zachary’s karate club, showed in Fig. 1, consists of the black community formed by a single node and the grey community formed by the rest of nodes. Although there is only one edge between these subsets, the grey community contains multiple subgroups of nodes with a high internal edge density, so it is possible to provide a different bipartition https://doi.org/10.1016/j.ejor.2024.02.033 Received 13 May 2023; Accepted 29 February 2024
European Journal of Operational Research 316 (2024) 519–538 520 D. Ponce et al. Fig. 1. Zachary’s karate club minimum cut. Fig. 2. Zachary’s karate club minimum normalized cut. of the Zachary’s karate club network with better community structure properties than the minimum cut bipartition. In order to solve this misconception, Shi and Malik (2000) introduced the normalized cut function aiming at solving some issues concerning the interpretability of the minimum cut problem. Normalized cuts compare the number of edges between different communities with the sum of the adjacency degree of each community, so the goal is to minimize the external edge density between nodes from different communities. For the basic case, given a graph 𝐺= (𝑉 , 𝐸)with vertex set 𝑉= {1,…, 𝑛}and edge set 𝐸, and a bipartition of the nodes set 𝐴, 𝐵 ⊂ 𝑉 , the normalized cut generated by {𝐴, 𝐵}is: ∑𝑖∈𝐴 𝑗∈𝐵 𝑤𝑖𝑗 ∑𝑖∈𝐴 𝑗∈𝑉 𝑤𝑖𝑗 +∑𝑖∈𝐴 𝑗∈𝐵 𝑤𝑖𝑗 ∑𝑖∈𝐵 𝑗∈𝑉 𝑤𝑖𝑗 ,(1) where 𝑤𝑖𝑗 is the number of edges between nodes 𝑖and 𝑗. We can define a function 𝑛𝑐𝑢𝑡(𝑆)for any subset 𝑆 ⊆ 𝑉 as: 𝑛𝑐𝑢𝑡(𝑆) = ∑𝑖∈𝑆 𝑗∉𝑆 𝑤𝑖𝑗 ∑𝑖∈𝑆 𝑗∈𝑉 𝑤𝑖𝑗 , so expression (1) is equal to 𝑛𝑐𝑢𝑡(𝐴) + 𝑛𝑐𝑢𝑡(𝐵). Thanks to 𝑤𝑖𝑗 definition, the normalized cut function can be applied to multigraphs where multiple edges are allowed. Example 1 (Cont.).Fig. 2 depicts the minimum normalized cut for the Zachary’s karate club graph. The normalized cut function (1) can be extended to weighted graphs whose weights are non-negative, finite and symmetric. These weights represent a similarity measure between nodes, where 𝑤𝑖𝑗 is the sum of the weights of the edges connecting 𝑖and 𝑗. As consequence, they are non-negative, finite and symmetric, i.e., 𝑤𝑖𝑗 =𝑤𝑗𝑖. Also, it can be assumed that 𝑤𝑖𝑖 = 0,∀𝑖∈𝑉, because generally the similarity of a node with itself is not considered. Nevertheless, for the solution methods introduced in this document, it would not be a problem. In addition, it is required that ∑𝑖∈𝐴 𝑗∈𝑉 𝑤𝑖𝑗 >0and ∑𝑖∈𝐵 𝑗∈𝑉 𝑤𝑖𝑗 >0for (1) to make sense. Normally, any graph that we want to partition must fulfill, without loss of generality, that ∑𝑗∈𝑉𝑤𝑖𝑗 >0for any 𝑖∈𝑉, because if this hypothesis is not fulfilled each node with an adjacent degree equal to zero can be considered as a community formed by itself, and therefore it would be removed from the graph to be considered. The previous normalized cut definition can be extended to the 𝑘way normalized cut, that partitions 𝑉in 𝑘non-empty subsets. Given a 𝑘-partition of the node set 𝐴1,…, 𝐴𝑘⊆ 𝑉 , the resulting extension of (1) to 𝑘-partitions is: 𝑘 ∑ 𝑠=1 ∑𝑖∈𝐴𝑠 𝑗∉𝐴𝑠 𝑤𝑖𝑗 ∑𝑖∈𝐴𝑠 𝑗∈𝑉 𝑤𝑖𝑗 = 𝑘 ∑ 𝑠=1 𝑛𝑐𝑢𝑡(𝐴𝑠).(2) The 𝑘-way normalized cut (2) is equivalent to sum the external edge density of each community of the 𝑘-partition. Shi and Malik (2000) propose to minimize the 2-way normalized cut function which depends on the characteristic vector of 2-partitions. These vectors are binary valued, thus the variables must be restricted to the binary domain to obtain the exact minimum normalized cut. Nonetheless, the proposed method, based on minimizing the Rayleigh quotient, finds the solution of the continuous relaxation of the problem as the eigenvector associated to the second smallest eigenvalue of the normalized Laplacian matrix. Thus, the solution is an approximation to the actual value. The normalized Laplacian matrix has been analyzed in multiple papers and books, such as Cheeger (1971), Donath and Hoffman (1973), Fiedler (1975), or Pothen et al. (1990), due to the interesting graph properties that it encodes. The paper by Shi and Malik (2000) proposed to apply 2-way normalized cuts recursively to obtain a final 𝑘-partition community structure. Obviously, this approach is clearly a heuristic since the best 𝑘-partition might not be found applying iteratively the 2-partition algorithm to previous elements of a partition. This is an important drawback of the eigenvector method that to date has not been solved. In the same paper, another heuristic algorithm is also presented to approximate the minimum 𝑘-way normalized cut based on greedy search and other spectral properties of the normalized Laplacian matrix. For example, they used eigenvectors associated to the smallest eigenvalues. In all cases, the optimality of the resulting partitions cannot be ensured. In other papers that revise the normalized cuts, as Barrientos and Madrid (2011), Cai et al. (2006), Hansen et al. (2012), He and Zhang (2016), Ojeda-Ruiz and Lee (2020), Tepper et al. (2011), Wang et al. (2022), Xu et al. (2009), Yang et al. (2022), or Zhong and Pun (2022), alternative algorithms have been developed aiming at obtaining community structures with a high internal edge density. These algorithms are based on the idea of relaxing the problem. However, no attempt has been made to obtain the exact minimum in any of the cases, and the idea of relaxing the problem is taken for granted. Our proposal is different from previous approaches in the literature. Rather than solving a relaxed version of the problem, we develop several reformulations using mathematical programming tools to exactly obtain the minimum 𝑘-way normalized cut of a graph. This is a challenging task since the normalized cut function is non-linear because it is defined by the ratio of two linear functions of integer variables. In spite of that, we are able to reformulate the minimum normalized cut problem as a Mixed-Integer Linear Programming (MILP) for which there are a variety of commercial solvers that can obtain the exact optimum. We note in passing that there are no previous exact approaches to solve this problem in the literature. Obviously, the problem is -complete (see Shi & Malik,2000), and trying to get its optimal solution makes the complexity of the methods to grow but still they are easily applicable for medium-sized graphs, although their running times are higher than other heuristic approaches, such as the spectral algorithms. The application of mathematical programming tools in community detection algorithms is not new, and one can find attempts to be used with several well-recognized measures. Among them, we cite the
European Journal of Operational Research 316 (2024) 519–538 521 D. Ponce et al. results in Costa (2015), Costa et al. (2017), Sukeda et al. (2023), where the modularity density measure is reformulated as a MILP yielding an improvement of the computation time for the exact maximization modularity density problem. We also cite the exact models to detect overlapping communities in networks proposed in Benati et al. (2022b, 2023). The contribution of this paper are: 1. To provide the first valid MILP formulation for the 𝑘-way normalized cut problem in graphs. 2. To develop a branch-and-price algorithm to solve the abovementioned problem. 3. To analyze extensive computational experiments that assess the usefulness of our methods to find the minimum 𝑘-way normalized cut. 4. To study a case study on segmentation of images showing the applicability of our algorithms to detect community structures in graphs. The rest of the paper is organized as follows: in Section 2, we formally define the minimum 𝑘-way normalized cut problem and propose some MILP formulations. In Section 3, a branch-and-price algorithm is defined to break the inherent symmetry of clustering problems. Next, in Section 4, the best proposed formulations and the branch-andprice algorithm are compared in an extensive computational analysis. Furthermore, the methods are applied to get segmentations of real images. Finally, the document ends with some conclusions in Section 5. 2. Some MILP formulations for the minimum 𝒌-way normalized cut problem In this section, we show how to formulate the normalized cut problem as a MILP. In particular, we focus on the 𝑘-way normalized cut because the original normalized cut (1) is a particular case of expression (2). As it happens in others community detection problems, there are several ways to formulate the normalized cut problem. For this reason we will present some alternative formulations to then compare their efficiency and their solution times. Due to the fact that the number of communities is fixed at 𝑘, many of our formulations are based on the ones proposed in Ales and Knippel (2020) for the 𝑘-partitioning problem. Given a weighted graph 𝐺= (𝑉 , 𝐸, 𝑤)with node set 𝑉= {1,…, 𝑛}, edge set 𝐸and weight set {𝑤𝑖𝑗 }(𝑖,𝑗)∈𝐸that are non-negative, finite and symmetric, we can redefine the set of weights between any nodes 𝑖, 𝑗 ∈𝑉as follows: 𝑤𝑖𝑗 =⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ 𝑤𝑖𝑗 ,if (𝑖, 𝑗) ∈ 𝐸, 𝑤𝑗𝑖,if (𝑗, 𝑖) ∈ 𝐸, 0,otherwise, so that, the minimum 𝑘-way normalized cut problem (MkWNCP) is defined as finding a 𝑘-partition {𝐴1,…, 𝐴𝑘}of 𝑉such that the objective function (2) is minimized. In order to formulate the problem as a MILP, we will follow a two steps procedure: (1) modeling the 𝑘-partition domain as a set of linear constraints and (2) modeling the non-linear expression (2) to be minimized as a MILP. Thus, we use variations of the formulations from Ales and Knippel (2020) to model the 𝑘-partition and then we will propose new formulations for the 𝑘-way normalized cut (2). First, we present two formulations for the 𝑘-partition, based on the assignment of each node to its community. Because of that, this type of formulations is called node-cluster. In both of them, we assume that the index of a community is always greater than or equal to the indices of the nodes that belong to it, reducing the number of variables. In addition, some conditions are imposed on the formulations to remove symmetric solutions representing the same 𝑘-partition. In the first formulation, it is considered that there exist 𝑘different non-empty communities whose indices belong to the set {𝑛−𝑘+1,…, 𝑛} and, as we noted above, a node cannot belong to a community if its index is greater than the community index. However, this does not remove all the symmetry of this formulation. Based on a formulation proposed in Carrizosa et al. (2020), we impose that the 𝑘communities 𝐴𝑛−𝑘+1,…, 𝐴𝑛are ordered by the maximum node index of each community, i.e., given two communities 𝑠and 𝑠′if 𝑛−𝑘+ 1 ≤𝑠<𝑠′≤𝑛, then max𝑖∈𝐴𝑠𝑖 < max𝑖∈𝐴𝑠′𝑖. This means that nodes are assigned to their communities from node 𝑛to node 1and communities are also filled in descending order from community 𝑛to 𝑛−𝑘+ 1. Due to the above hypothesis, the final 𝑘-partition is defined by the variables 𝑧𝑖𝑠 that are equivalent to: 𝑧𝑖𝑠 =⎧ ⎪ ⎨ ⎪ ⎩ 1,if node 𝑖is assigned to community 𝑠, 0,otherwise, for any 𝑠=𝑛−𝑘+ 1,…, 𝑛 and 𝑖= 1,…, 𝑠, and the first node-cluster formulation considered for the 𝑘-partition is (𝐹𝑎 𝑛𝑐 )∶ 𝑛 ∑ 𝑠=max{𝑖,𝑛−𝑘+1} 𝑧𝑖𝑠 = 1,∀𝑖= 1,…, 𝑛, (3) 𝑠 ∑ 𝑖=1 𝑧𝑖𝑠 ≥1,∀𝑠=𝑛−𝑘+ 1,…, 𝑛, (4) 𝑛 ∑ 𝑙=𝑠 𝑧𝑖𝑙 + 𝑠 ∑ 𝑗=𝑖+1 𝑧𝑗𝑠 ≥1,∀𝑠=𝑛−𝑘+ 2,…, 𝑛 − 1, 𝑖 = 2,…, 𝑠 − 1, (5) 𝑧𝑖𝑠 ∈ {0,1},∀𝑠=𝑛−𝑘+ 1,…, 𝑛, 𝑖 = 1,…, 𝑠. (6) Inequalities (3)–(4) impose that each node belongs to one and only one community and the number of non-empty communities is equal to 𝑘, respectively. Constraints (5) ensures the communities ordering that we have explained above. More specifically, if the highest node index of a community 𝑠is lower than the highest one of a community 𝑠′, then 𝑠<𝑠′. This condition needs not to be imposed when 𝑠=𝑛, 𝑠=𝑛−𝑘+ 1,𝑖=𝑠, or 𝑖= 1 because, in those cases, inequalities (3)–(4) ensure the condition. Finally, variables 𝑧𝑖𝑠 are defined as binary by (6). In general, removing all symmetric solutions leads to an improvement of the solution time. While in the first formulation it is imposed that the number of communities is equal to 𝑘and all of them are non-empty, in the second one it is considered that there exist 𝑛communities and exactly 𝑘of them are non-empty. Then, each non-empty community index is represented by the highest index of the nodes that belong to the community. This assumptions can be found in other papers, as Ales and Knippel (2020) or Benati et al. (2017), and allow us to remove all symmetric solutions of the problem that represent the same partition. Due to the above hypothesis, we can use the same 𝑧-variables but changing the indices domain. In this case, 𝑧𝑖𝑠 is defined for any 𝑠= 1,…, 𝑛, and 𝑖= 1,…, 𝑠. So, the second node-cluster formulation consists of: (𝐹𝑏 𝑛𝑐 )∶ 𝑛 ∑ 𝑠=𝑖 𝑧𝑖𝑠 = 1,∀𝑖= 1,…, 𝑛, (7) 𝑧𝑖𝑠 ≤𝑧𝑠𝑠,∀𝑖= 1,…, 𝑠 − 1, 𝑠 = 2,…, 𝑛 − 1,(8) 𝑛 ∑ 𝑠=1 𝑧𝑠𝑠 =𝑘, (9) 𝑧𝑖𝑠 ∈ {0,1},∀𝑠= 1,…, 𝑛, 𝑖 = 1,…, 𝑠. (10) Inequalities (7)–(9) impose that each node belongs to one and only one community, the highest node index of a non-empty community is equal to the community index and the number of non-empty communities is equal to 𝑘, respectively. Variables 𝑧𝑖𝑠 are defined as binary by (10).
European Journal of Operational Research 316 (2024) 519–538 522 D. Ponce et al. The parameter 𝑘has a clear influence in both formulations, especially in the first one, where a large value of 𝑘causes a considerable increase in the number of variables. So, when 𝑘is large enough the difference between the number of variables in both formulations is not significant. However, the performance regarding the bounds provided by their linear relaxations is empirically better in the second one as can be seen in Section 4. This explains that our formulations based on (𝐹𝑎 𝑛𝑐 ) report better performance for small values of 𝑘whereas the ones based on (𝐹𝑏 𝑛𝑐 ) are the best ones for larger values of 𝑘, as shown in our computational experience. Once we are able to model 𝑘-partitions of any graph thanks to (𝐹𝑎 𝑛𝑐 ) and (𝐹𝑏 𝑛𝑐 ), we can combine them with other elements to model the non-linear terms of the 𝑘-way normalized cut (2). We propose two types of 𝑘-way normalized cut representation, presented in the following subsections, and for each type we will build the corresponding formulations. 2.1. Three-index formulations In this subsection, we focus on the first type of formulation that uses variables with three indices. In particular, we develop two formulations of this type based on the domains (𝐹𝑎 𝑛𝑐 ) and (𝐹𝑏 𝑛𝑐 ), respectively. Under the same hypothesis as for the domain (𝐹𝑎 𝑛𝑐 ), the next two families of continuous variables: 𝛼𝑖𝑠 ≥0for any 𝑠=𝑛−𝑘+ 1,…, 𝑛 and 𝑖= 1,…, 𝑠; and 𝛽𝑖𝑗𝑠 ≥0for any 𝑠=𝑛−𝑘+ 1,…, 𝑛 and 𝑖, 𝑗 = 1,…, 𝑠, such that 𝑖 < 𝑗, are introduced to model the non-linear terms. They must fulfill the following conditions: 𝛼𝑖𝑠 𝑠 ∑ 𝑖′=1 𝑧𝑖′𝑠 𝑛 ∑ 𝑗′=1 𝑤𝑖′𝑗′=𝑧𝑖𝑠,∀𝑠=𝑛−𝑘+ 1,…, 𝑛, 𝑖 = 1,…, 𝑠, (11) 𝛽𝑖𝑗𝑠 𝑠 ∑ 𝑖′=1 𝑧𝑖′𝑠 𝑛 ∑ 𝑗′=1 𝑤𝑖′𝑗′=𝑧𝑖𝑠𝑧𝑗𝑠,∀𝑠=𝑛−𝑘+ 1,…, 𝑛, 𝑖, 𝑗 = 1,…, 𝑠, 𝑖 < 𝑗. (12) So, 𝛼𝑖𝑠 is equal to the inverse of the sum of node degrees from community 𝑠if node 𝑖belongs to 𝐴𝑠, and zero otherwise. Likewise, 𝛽𝑖𝑗𝑠 is equal to the inverse of the sum of node degrees from community 𝑠if nodes 𝑖and 𝑗belong to 𝐴𝑠. Otherwise, it is equal to zero. Assuming that {𝐴𝑛−𝑘+1,…, 𝐴𝑛}is a 𝑘-partition, these variables allow us to express the 𝑘-way normalized cut function (2) as: 𝑛 ∑ 𝑠=𝑛−𝑘+1 (𝑠−1 ∑ 𝑖=1 𝑠 ∑ 𝑗=𝑖+1 𝑤𝑖𝑗 (𝛼𝑖𝑠 +𝛼𝑗𝑠 − 2𝛽𝑖𝑗𝑠) + 𝑠 ∑ 𝑖=1 𝑛 ∑ 𝑗=𝑠+1 𝑤𝑖𝑗 𝛼𝑖𝑠),(13) because 𝛼𝑖𝑠 +𝛼𝑗𝑠 − 2𝛽𝑖𝑗𝑠 is equal to zero if 𝑖, 𝑗 ∈𝐴𝑠or 𝑖, 𝑗 ∉𝐴𝑠and it is equal to the inverse of the sum of node degrees from community 𝑠if 𝑖∈𝐴𝑠and 𝑗∉𝐴𝑠or 𝑗∈𝐴𝑠and 𝑖∉𝐴𝑠. In addition, if 𝑖≤𝑠<𝑗then 𝑗∉𝐴𝑠, so 𝛼𝑖𝑠 is equal to zero if 𝑖∉𝐴𝑠and it is equal to the inverse of the sum of node degrees from community 𝑠if 𝑖∈𝐴𝑠. Another important question to be considered, in order to define the relationships between variables, is which are the upper bounds of these continuous variables. It is not difficult to see that the maximum value that 𝛼𝑖𝑠 can take is 1 ∑𝑛 𝑗′=1 𝑤𝑖𝑗′when 𝐴𝑠= {𝑖}. The maximum value that 𝛽𝑖𝑗𝑠 can take is 1 ∑𝑛 𝑗′=1(𝑤𝑖𝑗′+𝑤𝑗𝑗′)when 𝐴𝑠= {𝑖, 𝑗}. Based on all the above ingredients, we propose the following formulation that uses variables with three indices: (𝐹𝑎 3𝑖)min 𝑛 ∑ 𝑠=𝑛−𝑘+1(𝑠−1 ∑ 𝑖=1 𝑠 ∑ 𝑗=𝑖+1 𝑤𝑖𝑗 (𝛼𝑖𝑠 +𝛼𝑗𝑠 − 2𝛽𝑖𝑗𝑠) + 𝑠 ∑ 𝑖=1 𝑛 ∑ 𝑗=𝑠+1 𝑤𝑖𝑗 𝛼𝑖𝑠) 𝑠.𝑡. ∶(3)–(6), 𝛼𝑖𝑠 ≤𝑧𝑖𝑠 1 ∑𝑛 𝑗′=1 𝑤𝑖𝑗′ ,∀𝑠=𝑛−𝑘+ 1,…, 𝑛, ∀𝑖= 1,…, 𝑠, (14) 𝛽𝑖𝑗𝑠 ≤𝑧𝑖𝑠 1 ∑𝑛 𝑗′=1(𝑤𝑖𝑗′+𝑤𝑗𝑗′),∀𝑠=𝑛−𝑘+ 1,…, 𝑛, ∀𝑖, 𝑗 = 1,…, 𝑠, 𝑖 < 𝑗, (15) 𝛽𝑖𝑗𝑠 ≤𝑧𝑗𝑠 1 ∑𝑛 𝑗′=1(𝑤𝑖𝑗′+𝑤𝑗𝑗′),∀𝑠=𝑛−𝑘+ 1,…, 𝑛, ∀𝑖, 𝑗 = 1,…, 𝑠, 𝑖 < 𝑗, (16) 𝛼𝑖𝑠 ≤𝛼𝑗𝑠 + (1 − 𝑧𝑗𝑠)1 ∑𝑛 𝑗′=1 𝑤𝑖𝑗′ ,∀𝑠=𝑛−𝑘+ 1,…, 𝑛, 𝑖, 𝑗 = 1,…, 𝑠, 𝑖 ≠𝑗, (17) 𝛽𝑖𝑗𝑠 ≤𝛼𝑖𝑠,∀𝑠=𝑛−𝑘+ 1,…, 𝑛, ∀𝑖, 𝑗 = 1,…, 𝑠, 𝑖 < 𝑗, (18) 𝛼𝑖𝑠 ≤𝛽𝑖𝑗𝑠 + (1 − 𝑧𝑗𝑠)1 ∑𝑛 𝑗′=1 𝑤𝑖𝑗′ ,∀𝑠=𝑛−𝑘+ 1,…, 𝑛, ∀𝑖, 𝑗 = 1,…, 𝑠, 𝑖 < 𝑗, (19) 𝛽𝑖𝑗𝑠 ≤𝛼𝑗𝑠,∀𝑠=𝑛−𝑘+ 1,…, 𝑛, ∀𝑖, 𝑗 = 1,…, 𝑠, 𝑖 < 𝑗, (20) 𝛼𝑗𝑠 ≤𝛽𝑖𝑗𝑠 + (1 − 𝑧𝑖𝑠)1 ∑𝑛 𝑗′=1 𝑤𝑗𝑗′ ,∀𝑠=𝑛−𝑘+ 1,…, 𝑛, ∀𝑖, 𝑗 = 1,…, 𝑠, 𝑖 < 𝑗, (21) 𝑠 ∑ 𝑖=1 𝛼𝑖𝑠 𝑛 ∑ 𝑗′=1 𝑤𝑖𝑗′= 1,∀𝑠=𝑛−𝑘+ 1,…, 𝑛, (22) 𝛼𝑖𝑠 ≥0,∀𝑠=𝑛−𝑘+ 1,…, 𝑛, ∀𝑖= 1,…, 𝑠, (23) 𝛽𝑖𝑗𝑠 ≥0,∀𝑠=𝑛−𝑘+ 1,…, 𝑛, ∀𝑖, 𝑗 = 1,…, 𝑠, 𝑖 < 𝑗. (24) Constraints (14)–(16) impose that 𝛼𝑖𝑠 = 0 if node 𝑖does not belong to community 𝑠, and 𝛽𝑖𝑗𝑠 = 0 if node 𝑖or node 𝑗does not belong to community 𝑠, i.e., 𝑧𝑖𝑠 = 0 or 𝑧𝑗𝑠 = 0. The families of constraints (17)–(21) guarantee that 𝛽𝑖𝑗𝑠 =𝛼𝑖𝑠 =𝛼𝑗𝑠 if nodes 𝑖and 𝑗belong to community 𝑠, i.e., 𝑧𝑖𝑠 =𝑧𝑗𝑠 = 1. Combining the above inequalities with equality (22), we ensure that variables 𝛼𝑖𝑠 and 𝛽𝑖𝑗𝑠 fulfill conditions (11) and (12). Finally, variables 𝛼𝑖𝑠 and 𝛽𝑖𝑗𝑠 are defined as non-negative and continuous by (23)–(24). The second three-index formulation based on (𝐹𝑏 𝑛𝑐 )uses the same variables 𝛼𝑖𝑠 and 𝛽𝑖𝑗𝑠 but changing their indices domain, since there are 𝑛communities and a community 𝑠is non-empty if and only if node 𝑠belongs to it. In this case, 𝛼𝑖𝑠 is defined for any 𝑠= 1,…, 𝑛 and 𝑖= 1,…, 𝑠, and 𝛽𝑖𝑗𝑠 is defined for any 𝑠= 1,…, 𝑛 and 𝑖, 𝑗 = 1,…, 𝑠−1 such that 𝑖<𝑗. These variables continue to fulfill the previous conditions (11) and (12) but for their new index set. The expression of the 𝑘-way normalized cut function (2) also changes in this alternative formulation: 𝑛 ∑ 𝑠=1 (𝑠−1 ∑ 𝑖=1 𝑤𝑖𝑠(𝛼𝑠𝑠 −𝛼𝑖𝑠) + 𝑠−2 ∑ 𝑖=1 𝑠−1 ∑ 𝑗=𝑖+1 𝑤𝑖𝑗 (𝛼𝑖𝑠 +𝛼𝑗𝑠 − 2𝛽𝑖𝑗𝑠) + 𝑠 ∑ 𝑖=1 𝑛 ∑ 𝑗=𝑠+1 𝑤𝑖𝑗 𝛼𝑖𝑠). (25) In this case, given a community 𝑠and a node 𝑖, such that 𝑖 < 𝑠, the expression 𝛼𝑠𝑠 −𝛼𝑖𝑠 is equal to zero if community 𝑠is empty or 𝑖∈𝐴𝑠, and it is equal to the inverse of the sum of node degrees from community 𝑠if community 𝑠is non-empty and 𝑖∉𝐴𝑠. The remaining terms of the expression (25) have been described in the objective function (13) of the above formulation (𝐹𝑎 3𝑖). The upper bounds of the variables are different: 𝛼𝑠𝑠 is bounded by 1 ∑𝑛 𝑗′=1 𝑤𝑠𝑗′for 𝐴𝑠= {𝑠};𝛼𝑖𝑠 is bounded by 1 ∑𝑛 𝑗′=1(𝑤𝑖𝑗′+𝑤𝑠𝑗′)for 𝐴𝑠= {𝑖, 𝑠}, 𝑖 < 𝑠; and 𝛽𝑖𝑗𝑠 is bounded by 1 ∑𝑛 𝑗′=1(𝑤𝑖𝑗′+𝑤𝑗𝑗′+𝑤𝑠𝑗′)for 𝐴𝑠= {𝑖, 𝑗, 𝑠}, 𝑖 < 𝑗 < 𝑠. Then, the second three-index formulation is the following one: (𝐹𝑏 3𝑖)min 𝑛 ∑ 𝑠=1 (𝑠−1 ∑ 𝑖=1 𝑤𝑖𝑠(𝛼𝑠𝑠 −𝛼𝑖𝑠) + 𝑠−2 ∑ 𝑖=1 𝑠−1 ∑ 𝑗=𝑖+1 𝑤𝑖𝑗 (𝛼𝑖𝑠 +𝛼𝑗𝑠 − 2𝛽𝑖𝑗𝑠) + 𝑠 ∑ 𝑖=1 𝑛 ∑ 𝑗=𝑠+1 𝑤𝑖𝑗 𝛼𝑖𝑠) 𝑠.𝑡. ∶(7)–(10), 𝛼𝑠𝑠 ≤𝑧𝑠𝑠 1 ∑𝑛 𝑗′=1 𝑤𝑠𝑗′ ,∀𝑠= 1,…, 𝑛, (26) 𝛼𝑖𝑠 ≤𝑧𝑖𝑠 1 ∑𝑛 𝑗′=1(𝑤𝑖𝑗′+𝑤𝑠𝑗′),∀𝑖, 𝑠 = 1,…, 𝑛, 𝑖 < 𝑠, (27) 𝛽𝑖𝑗𝑠 ≤𝑧𝑖𝑠 1 ∑𝑛 𝑗′=1(𝑤𝑖𝑗′+𝑤𝑗𝑗′+𝑤𝑠𝑗′),∀𝑖, 𝑗, 𝑠 = 1,…, 𝑛, 𝑖 < 𝑗 < 𝑠, (28)
European Journal of Operational Research 316 (2024) 519–538 523 D. Ponce et al. 𝛽𝑖𝑗𝑠 ≤𝑧𝑗𝑠 1 ∑𝑛 𝑗′=1(𝑤𝑖𝑗′+𝑤𝑗𝑗′+𝑤𝑠𝑗′),∀𝑖, 𝑗, 𝑠 = 1,…, 𝑛, 𝑖 < 𝑗 < 𝑠, (29) 𝛼𝑖𝑠 ≤𝛼𝑠𝑠,∀𝑖, 𝑠 = 1,…, 𝑛, 𝑖 < 𝑠, (30) 𝛼𝑠𝑠 ≤𝛼𝑖𝑠 + (1 − 𝑧𝑖𝑠)1 ∑𝑛 𝑗′=1 𝑤𝑠𝑗′ ,∀𝑖, 𝑠 = 1,…, 𝑛, 𝑖 < 𝑠, (31) 𝛽𝑖𝑗𝑠 ≤𝛼𝑖𝑠,∀𝑖, 𝑗, 𝑠 = 1,…, 𝑛, 𝑖 < 𝑗 < 𝑠, (32) 𝛼𝑖𝑠 ≤𝛽𝑖𝑗𝑠 + (1 − 𝑧𝑗𝑠)1 ∑𝑛 𝑗′=1(𝑤𝑖𝑗′+𝑤𝑠𝑗′),∀𝑖, 𝑗, 𝑠 = 1,…, 𝑛, 𝑖 < 𝑗 < 𝑠, (33) 𝛽𝑖𝑗𝑠 ≤𝛼𝑗𝑠,∀𝑖, 𝑗, 𝑠 = 1,…, 𝑛, 𝑖 < 𝑗 < 𝑠, (34) 𝛼𝑗𝑠 ≤𝛽𝑖𝑗𝑠 + (1 − 𝑧𝑖𝑠)1 ∑𝑛 𝑗′=1(𝑤𝑗𝑗′+𝑤𝑠𝑗′),∀𝑖, 𝑗, 𝑠 = 1,…, 𝑛, 𝑖 < 𝑗 < 𝑠, (35) 𝑠 ∑ 𝑖=1 𝛼𝑖𝑠 𝑛 ∑ 𝑗′=1 𝑤𝑖𝑗′=𝑧𝑠𝑠,∀𝑠= 1,…, 𝑛, (36) 𝛼𝑖𝑠 ≥0,∀𝑖, 𝑠 = 1,…, 𝑛, 𝑖 ≤𝑠, (37) 𝛽𝑖𝑗𝑠 ≥0,∀𝑖, 𝑗, 𝑠 = 1,…, 𝑛, 𝑖 < 𝑗 < 𝑠. (38) Constraints (26)–(29) impose that 𝛼𝑠𝑠 is equal to zero if community 𝑠is empty, 𝛼𝑖𝑠 is equal to zero if node 𝑖does not belong to community 𝑠and 𝛽𝑖𝑗𝑠 is equal to zero if node 𝑖or 𝑗does not belong to community 𝑠. The families of constraints (30)–(35) guarantee that for each community 𝑠the non-null variables 𝛼𝑖𝑠 and 𝛽𝑖𝑗𝑠 are equal. Because of this and equality (36), we ensure that, for each non-empty community 𝑠, each non-null variable 𝛼𝑖𝑠 or 𝛽𝑖𝑗𝑠 is equal to the inverse of the sum of node degrees from community 𝑠. Finally, variables 𝛼𝑖𝑠 and 𝛽𝑖𝑗𝑠 are defined as non-negative and continuous in (37)–(38). Despite having proposed two alternative exact MILP formulations for the 𝑘-way normalized cut, (𝐹𝑎 3𝑖) and (𝐹𝑏 3𝑖) can be simplified because, even if we remove some of their constraints, the optimal value remains the same. This is the case of inequalities (4),(19), and (21) of (𝐹𝑎 3𝑖). So, if (𝐹𝑎∗ 3𝑖) is the model that results from (𝐹𝑎 3𝑖) by removing constraints (4),(19) and (21), we can prove the following proposition. Proposition 1. Let 𝑓𝑎 3𝑖and 𝑓𝑎∗ 3𝑖be the optimal values of (𝐹𝑎 3𝑖) and (𝐹𝑎∗ 3𝑖), respectively. Then: 𝑓𝑎 3𝑖=𝑓𝑎∗ 3𝑖. Proof. First, (4) is obtained by applying inequality (14) over each 𝛼𝑖𝑠 of equality (22). In case that 𝑧𝑖𝑠 =𝑧𝑗𝑠 = 1,𝛽𝑖𝑗𝑠 could be nonnull and, constraints (17),(18) and (20) impose that 𝛽𝑖𝑗𝑠 ≤𝛼𝑖𝑠 =𝛼𝑗𝑠. Since constraints (19) and (21) only impose lower bounds for 𝛽𝑖𝑗𝑠 variables, we can find an optimal solution of (𝐹𝑎∗ 3𝑖) where 𝛽𝑖𝑗𝑠 =𝛼𝑖𝑠 = 𝛼𝑗𝑠 because 𝛽𝑖𝑗𝑠 is always multiplied by a non-positive value in the objective function and the model forces to increase the value of this kind of variables as much as possible, which means that lower bounds have no influence. □ We can follow a similar argument with the alternative formulation (𝐹𝑏 3𝑖). So, considering that (𝐹𝑏∗ 3𝑖) is the model that results from (𝐹𝑏 3𝑖) by removing constraints (26),(33), and (35) we can state the following proposition. Proposition 2. Let 𝑓𝑏 3𝑖and 𝑓𝑏∗ 3𝑖be the optimal values of (𝐹𝑏 3𝑖) and (𝐹𝑏∗ 3𝑖), respectively. Then: 𝑓𝑏 3𝑖=𝑓𝑏∗ 3𝑖. Proof. As in the proof of Proposition 1, inequalities (33) and (35) impose lower bounds for 𝛽-variables, a family of variables that appear in the objective function with negative sign, thus they will assume values as larger as possible making the lower bounds to be useless. Hence, these constraints have no influence on the optimal solution of model (𝐹𝑏 3𝑖) and they become redundant. Finally, if an optimal solution satisfies equality (36), then this condition implies directly inequality (26), since: 𝑧𝑠𝑠 = 𝑠 ∑ 𝑖=1 𝛼𝑖𝑠 𝑛 ∑ 𝑗′=1 𝑤𝑖𝑗′≥𝛼𝑠𝑠 𝑛 ∑ 𝑗′=1 𝑤𝑠𝑗′⇒𝛼𝑠𝑠 ≤𝑧𝑠𝑠 1 ∑𝑛 𝑗′=1 𝑤𝑠𝑗′ .□ There are some other constraints that are also redundant in the MILP formulation as, for instance, (15),(16),(28), and (29). However, we keep them in formulations (𝐹𝑎∗ 3𝑖) and (𝐹𝑏∗ 3𝑖) because removing them worsens the lower bound of the continuous relaxation, which affects the solution time of the actual problem with discrete variables. Also, we have analyzed the specific case of the problem when 𝑘= 2 because for that important case we can shed more light into the structure of the solutions. In particular, we can state the following proposition. Proposition 3. Given a weighted graph 𝐺= (𝑉 , 𝐸, 𝑤)with non-negative symmetric weights and assuming that {𝐴, 𝐵}is the minimum 2-way cut of 𝐺, 𝐴cannot strictly contain a subset with a normalized cut value strictly lower than 𝑛𝑐𝑢𝑡(𝐴)and 𝐵cannot strictly contain a subset with a normalized cut value strictly lower than 𝑛𝑐𝑢𝑡(𝐵). Proof. We will just prove the proposition for 𝐴because for 𝐵is analogous. Assuming that there is a subset 𝐴′such that 𝐴′⊂ 𝐴 and 𝑛𝑐𝑢𝑡(𝐴′)< 𝑛𝑐𝑢𝑡(𝐴), it follows that: ∑𝑖∈𝐴′ 𝑗∈𝐴⧵𝐴′ 𝑤𝑖𝑗 +∑𝑖∈𝐴′ 𝑗∈𝐵 𝑤𝑖𝑗 ∑𝑖∈𝐴′ 𝑗∈𝑉 𝑤𝑖𝑗 =𝑛𝑐𝑢𝑡(𝐴′)< 𝑛𝑐𝑢𝑡(𝐴) =∑𝑖∈𝐴′ 𝑗∈𝐵 𝑤𝑖𝑗 +∑𝑖∈𝐴⧵𝐴′ 𝑗∈𝐵 𝑤𝑖𝑗 ∑𝑖∈𝐴 𝑗∈𝑉 𝑤𝑖𝑗 <∑𝑖∈𝐴′ 𝑗∈𝐵 𝑤𝑖𝑗 +∑𝑖∈𝐴⧵𝐴′ 𝑗∈𝐵 𝑤𝑖𝑗 ∑𝑖∈𝐴′ 𝑗∈𝑉 𝑤𝑖𝑗 ⇒ ∑ 𝑖∈𝐴′ 𝑗∈𝐴⧵𝐴′ 𝑤𝑖𝑗 +∑ 𝑖∈𝐴′ 𝑗∈𝐵 𝑤𝑖𝑗 <∑ 𝑖∈𝐴′ 𝑗∈𝐵 𝑤𝑖𝑗 +∑ 𝑖∈𝐴⧵𝐴′ 𝑗∈𝐵 𝑤𝑖𝑗 ⇒∑ 𝑖∈𝐴′ 𝑗∈𝐴⧵𝐴′ 𝑤𝑖𝑗 <∑ 𝑖∈𝐴⧵𝐴′ 𝑗∈𝐵 𝑤𝑖𝑗 . Now, defining a new subset 𝐵′=𝐵∪ (𝐴⧵𝐴′)we can observe the following relationship between the normalized cuts of 𝐵and 𝐵′: 𝑛𝑐𝑢𝑡(𝐵′) = ∑𝑖∈𝐵 𝑗∈𝐴′𝑤𝑖𝑗 +∑𝑖∈𝐴⧵𝐴′ 𝑗∈𝐴′ 𝑤𝑖𝑗 ∑𝑖∈𝐵′ 𝑗∈𝑉 𝑤𝑖𝑗 <∑𝑖∈𝐵 𝑗∈𝐴′𝑤𝑖𝑗 +∑𝑖∈𝐵 𝑗∈𝐴⧵𝐴′𝑤𝑖𝑗 ∑𝑖∈𝐵′ 𝑗∈𝑉 𝑤𝑖𝑗 <∑𝑖∈𝐵 𝑗∈𝐴 𝑤𝑖𝑗 ∑𝑖∈𝐵 𝑗∈𝑉 𝑤𝑖𝑗 =𝑛𝑐𝑢𝑡(𝐵). Hence, the 2-way normalized cut obtained by the bipartition {𝐴′, 𝐵′}is strictly lower than the one obtained by {𝐴, 𝐵}which is a contradiction and we can ensure that 𝐴cannot strictly contain a subset with a normalized cut value strictly lower than 𝑛𝑐𝑢𝑡(𝐴).□ We will use this last observation to prove that we can remove inequalities (17) from (𝐹𝑎∗ 3𝑖) and (31) from (𝐹𝑏∗ 3𝑖) when 𝑘= 2, which imply that for each community 𝑠and, nodes 𝑖and 𝑗, such that 𝑖, 𝑗 ∈𝐴𝑠, variables 𝛼𝑖𝑠 and 𝛼𝑗𝑠 are equal, but first we have to prove the following proposition. Proposition 4. Given a weighted graph 𝐺= (𝑉 , 𝐸, 𝑤)with non-negative symmetric weights, a non-empty subset 𝑆 ⊆ 𝑉 , a vector of positive values {𝑉𝑖}𝑖∈𝑆and a constant 𝑎≥0, the optimal value of: (𝑀1)max ∑ 𝑖,𝑗∈𝑆 𝑤𝑖𝑗 min{𝑎𝑖, 𝑎𝑗}(39) ∑ 𝑖∈𝑆 𝑎𝑖𝑉𝑖=𝑎, (40) 𝑎𝑖≥0,∀𝑖∈𝑆, (41) is equal to: 𝑎⋅max {∑𝑖,𝑗∈𝑆′𝑤𝑖𝑗 ∑𝑖∈𝑆′𝑉𝑖 ∶ ∅ ≠𝑆′⊆ 𝑆},(42)
European Journal of Operational Research 316 (2024) 519–538 530 D. Ponce et al. 4.1. Comparing formulations performance This subsection assesses the performance of each one of the exact formulations proposed in Section 2: (𝐹𝑎 3𝑖), (𝐹𝑏 3𝑖), (𝐹𝑎∗ 3𝑖), (𝐹𝑏∗ 3𝑖), (𝐹𝑎+ 3𝑖), (𝐹𝑏+ 3𝑖), (𝐹𝑎∗ 2𝑖), and (𝐹𝑏∗ 2𝑖). Based on Ales and Knippel (2020), we have generated 200 undirected weighted graphs without loops which weights are randomly drawn using a uniform distribution in the interval [0,500]: 100 graphs of 10 nodes and 100 graphs of 15 nodes. We have solved these instances for different values of the parameter 𝑘∈ {2,3,4,5,6,7, 8,9}. We report average results of all the instances for different combinations of 𝑛and 𝑘. Due to the large number of considered random instances we set a time limit of 30 minutes on the solution time for solving each of them. This implies that, in many cases, we cannot certify optimality for the solutions provided by some of the formulations. The results of the above experiments are reported in Tables 1 and 2, and they correspond to the average for each combination of parameters 𝑛and 𝑘. We have reported the running times in seconds (Time (s)), the gap (100(upper bound −lower bound)∕upper bound) at termination (GAP (%)), the number of nodes explored in the branchand-bound (B&B) tree obtained by the proposed formulations (Nodes), the percentage of instances solved to optimality (Opt (%)) and the percentage of instances in which the root node was solved (Dual Bound (%)). It is important to note that, in case there is no other lower bound than zero when the linear relaxation could not be solved at root node, it is reported a gap equal to 100%. First of all, in Table 1, we have compared formulations (𝐹𝑎∗ 3𝑖) and (𝐹𝑏∗ 3𝑖) with their original versions (𝐹𝑎 3𝑖) and (𝐹𝑏 3𝑖) for 𝑘= 2, to show the improvement obtained by Propositions 1–5. Although the difference is not very significant, formulation (𝐹𝑎∗ 3𝑖) performs better than (𝐹𝑎 3𝑖) and (𝐹𝑏∗ 3𝑖) better than (𝐹𝑏 3𝑖) in terms of time, for every size 𝑛. So, we can conclude that the improvement is clear and we will only consider the reduced formulations (𝐹𝑎∗ 3𝑖) and (𝐹𝑏∗ 3𝑖) in the following studies, since they are simpler than their original versions. Looking at the results of Table 2, there is a clear evidence that the best performing formulations are (𝐹𝑎+ 3𝑖) and (𝐹𝑏+ 3𝑖) thanks to the improvement brought about by equalities (44) and (45). The reader can observe in that table that all our formulations are able to solve to optimality all instances with 𝑛= 10 in less than 1800 seconds, and the best CPU-times are obtained by (𝐹𝑎+ 3𝑖) and (𝐹𝑏+ 3𝑖). On the contrary, for 𝑛= 15, there is a number of instances where optimality of the solutions found cannot be certified using formulations (𝐹𝑎∗ 3𝑖), (𝐹𝑏∗ 3𝑖), (𝐹𝑎∗ 2𝑖), or (𝐹𝑏∗ 2𝑖), whereas using (𝐹𝑎+ 3𝑖) or (𝐹𝑏+ 3𝑖) all instances are solved up to optimality. More specifically, according with our experiments, the worst performance is obtained with formulations (𝐹𝑎∗ 2𝑖) and (𝐹𝑏∗ 2𝑖) and the reason may be the low quality of the bounds provided by their continuous relaxations, even though the number of variables is much smaller. In general, the number of branching nodes of (𝐹𝑎+ 3𝑖) and (𝐹𝑏+ 3𝑖) is lower than the one obtained by the remaining formulations. (𝐹𝑏+ 3𝑖) explores a slightly fewer number of nodes because formulation (𝐹𝑏 𝑛𝑐 ) enjoys better continuous relaxation bounds. On the contrary, (𝐹𝑎∗ 2𝑖) and (𝐹𝑏∗ 2𝑖) need to explore more nodes by the worse quality of their continuous relaxation bounds. Regarding the value of the parameter 𝑘, we observe that formulations based on (𝐹𝑎 𝑛𝑐 ) perform better for small values of 𝑘since they are the ones with the smaller number of variables even though their bounds are worse. However, if 𝑘increases this effect is mitigated, and the influence of the bounds compensates the dimension of the formulations so that those based on (𝐹𝑏 𝑛𝑐 ) become the best. The performance of (𝐹𝑎∗ 2𝑖) and (𝐹𝑏∗ 2𝑖) does not improve for bigger instances (𝑛= 15). This behavior leads us not to consider these formulations in the following. Finally, in order to analyze in more detail the linear relaxations of our formulations, we report the average gap at the root node Table 1 Computational analysis for 𝑘= 2. Formulation 𝑛 10 15 20 𝐹𝑎 3𝑖Time (s) 0.16 2.85 79.89 GAP (%) 0.00 0.00 0.00 Nodes 297 5107 109 808 𝐹𝑏 3𝑖Time (s) 0.46 7.11 181.39 GAP (%) 0.00 0.00 0.00 Nodes 183 5017 100 449 𝐹𝑎∗ 3𝑖Time (s) 0.14 2.41 77.00 GAP (%) 0.00 0.00 0.00 Nodes 300 5531 114 880 𝐹𝑏∗ 3𝑖Time (s) 0.46 6.45 181.00 GAP (%) 0.00 0.00 0.00 Nodes 185 4996 101 340 (100(optimal value −linear relaxation)∕optimal value) in Table 3. As discussed in the previous sections, the formulations with best linear relaxations are (𝐹𝑎+ 3𝑖) and (𝐹𝑏+ 3𝑖), then (𝐹𝑎∗ 3𝑖) and (𝐹𝑏∗ 3𝑖), and the worst linear relaxations are obtained by (𝐹𝑎∗ 2𝑖) and (𝐹𝑏∗ 2𝑖). Moreover, the formulations based on (𝐹𝑏 𝑛𝑐 ) obtain better linear relaxation results compared to the ones based on (𝐹𝑎 𝑛𝑐 ). This fact can be observed as well in the number of explored nodes of each formulation in Table 2. 4.2. Random instances Also based on Ales and Knippel (2020), we have built 200 undirected weighted graphs without loops whose weights are randomly generated using an uniform distribution in the interval [0,500]:50 graphs for each size 𝑛∈ {20,30,40,50}. For each family of random graphs, we have considered values of parameter 𝑘∈ {2,4,6,8}. The formulations used in the experiments are those that performed better in the preliminary computational experiments in Section 4.1: (𝐹𝑎∗ 3𝑖), (𝐹𝑏∗ 3𝑖), (𝐹𝑎+ 3𝑖), and (𝐹𝑏+ 3𝑖). These formulations are compared with the branch-and-price algorithm. Due to the big amount of random instances we test, we consider a relatively small maximal computational time of 30 minutes, so most of the instances with more than 30 nodes are not solved to optimality. The fact that all the instances are complete graph with random weights in the interval [0,500] provokes that these random graphs are more complex than other real networks that we can find in literature. For the implementation of the B&P procedure we have applied a multiple pricing strategy, adding several new columns in each iteration of the column generation routine. Specifically, for 𝑛= 20, we add up to a maximum number of 10 feasible columns with negative reduced cost, for 𝑛= 30, a maximum number of 30, for 𝑛= 40, a maximum number of 40 and, for 𝑛= 50, a maximum number 50. The average gap at the root node values are reported in Table 4. In this case, we have just considered instances that are solved at optimality, i.e., when 𝑛∈ {20,30}. The results in Table 4 point out that the best method in terms of linear relaxation is our B&P procedure. Better linear relaxation value will lead eventually to fewer nodes explored. The results of the above experiments are reported in Table 5. As before, they correspond to the average running time, average gap at termination, average number nodes explored in the B&B tree, the percentage of instances solved to optimality and the percentage of instances in which the root node was solved after solving 50 instances for each combination of parameters 𝑛and 𝑘. Looking at Table 5, we can conclude that for the smallest value of 𝑘 (𝑘= 2) the best results are obtained by formulation (𝐹𝑎+ 3𝑖), that may be influenced by its reduced number of variables. When 𝑘increases, the best results are obtained with the B&P procedure. Formulations based on (𝐹𝑎 𝑛𝑐 ) turn out to be the worst ones. Furthermore, we can conclude that the B&P procedure works really well since it is the only method able to exactly solve all the instances
European Journal of Operational Research 316 (2024) 519–538 531 D. Ponce et al. Table 2 Computational results of formulations in Section 2for small-sized instances. 𝑛 10 15 Formulation 𝑘 2 3 4 5 6 7 8 9 2 3 4 5 6 7 8 9 𝐹𝑎∗ 3𝑖Time (s) 0.14 1.15 2.91 3.33 2.39 1.31 0.71 0.10 2.41 163.89 1544.92 1800.00 1800.00 1800.00 1750.00 924.67 GAP (%) 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 6.70 15.10 13.00 9.40 4.20 0.00 Nodes 300 2171 3772 3330 1779 454 45 1 5531 272 634 1 474 439 1030 093 747 843 676565 855 309 512 951 Opt (%) 100 100 100 100 100 100 100 100 100 100 53 0 0 1 20 100 Dual Bound (%) 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 𝐹𝑏∗ 3𝑖Time (s) 0.46 1.52 1.36 1.02 0.72 0.51 0.24 0.03 6.45 227.00 837.61 898.00 579.90 309.36 166.90 91.33 GAP (%) 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.07 0.08 0.00 0.00 0.00 0.00 Nodes 185 1013 1016 701 403 124 14 1 4996 144135 397 796 403 172 277052 162 695 94260 47235 Opt (%) 100 100 100 100 100 100 100 100 100 100 99 98 100 100 100 100 Dual Bound (%) 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 𝐹𝑎+ 3𝑖Time (s) 0.04 0.09 0.08 0.06 0.05 0.03 0.03 0.04 0.31 1.54 1.92 1.72 1.18 0.77 0.54 0.41 GAP (%) 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 Nodes 4 5 2 1 1 1 1 1 69 165 127 56 19 3 2 1 Opt (%) 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 Dual Bound (%) 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 𝐹𝑏+ 3𝑖Time (s) 0.17 0.12 0.08 0.05 0.03 0.02 0.02 0.02 3.41 4.13 2.11 1.07 0.49 0.28 0.19 0.15 GAP (%) 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 Nodes 2 1 1 1 1 1 1 1 25 54 16 4 2 1 1 1 Opt (%) 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 Dual Bound (%) 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 𝐹𝑎∗ 2𝑖Time (s) 0.31 3.35 9.21 12.64 8.02 2.48 0.66 0.21 6.19 442.20 1800.00 1800.00 1800.00 1800.00 1800.00 1800.00 GAP (%) 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 42.90 59.20 59.10 51.70 41.80 28.00 Nodes 794 11 216 35 871 45 042 26 191 7777 513 2 8604 877 049 3 467 219 2 739 554 2 738 054 2 904 745 3 175 638 3 915 061 Opt (%) 100 100 100 100 100 100 100 100 100 100 0 0 0 0 0 0 Dual Bound (%) 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 𝐹𝑏∗ 2𝑖Time (s) 0.37 3.91 10.06 13.50 8.26 2.34 0.56 0.14 10.90 520.30 1800.00 1800.00 1800.00 1800.00 1800.00 1800.00 GAP (%) 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 58.90 61.20 64.80 43.30 40.50 14.90 Nodes 598 11 113 34 447 42 605 24 620 6924 345 1 13 161 822 065 2 898 647 2 534 913 2 744 306 2 469 526 3 316 437 3 866 377 Opt (%) 100 100 100 100 100 100 100 100 100 100 0 0 0 0 0 0 Dual Bound (%) 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 Table 3 Gap at the root node for small-sized instances. 𝑛Formulation 𝑘 2 3 456789 10 𝐹𝑎∗ 3𝑖71.29 73.67 62.04 45.17 33.31 21.90 9.90 0.02 𝐹𝑏∗ 3𝑖66.67 51.83 30.43 18.22 11.09 5.93 2.35 0.00 𝐹𝑎+ 3𝑖2.73 1.27 0.15 0.00 0.00 0.00 0.00 0.00 𝐹𝑏+ 3𝑖1.01 0.13 0.02 0.00 0.00 0.00 0.00 0.00 𝐹𝑎∗ 2𝑖100.00 100.00 95.88 75.89 71.15 89.02 94.68 10.12 𝐹𝑏∗ 2𝑖100.00 100.00 94.48 74.58 67.49 96.88 84.62 2.28 15 𝐹𝑎∗ 3𝑖81.36 59.45 62.46 61.66 48.16 33.88 25.47 20.19 𝐹𝑏∗ 3𝑖79.57 57.38 46.36 27.01 20.03 15.37 11.34 8.20 𝐹𝑎+ 3𝑖14.28 6.54 2.95 1.05 0.30 0.06 0.01 0.00 𝐹𝑏+ 3𝑖13.97 3.86 1.07 0.22 0.05 0.00 0.00 0.00 𝐹𝑎∗ 2𝑖100.00 89.68 94.20 95.99 97.14 97.78 98.14 97.50 𝐹𝑏∗ 2𝑖99.22 90.90 96.40 96.63 98.29 97.45 96.66 89.17 Table 4 Gap at the root node for medium-sized instances. 𝑛Formulation 𝑘 2468 20 𝐹𝑎∗ 3𝑖56.73 63.63 72.46 68.56 𝐹𝑏∗ 3𝑖60.41 57.99 46.57 29.05 𝐹𝑎+ 3𝑖21.23 5.78 1.72 0.46 𝐹𝑏+ 3𝑖23.09 3.61 0.72 0.11 B&P 0.00 0.02 0.02 0.00 30 𝐹𝑎∗ 3𝑖75.44 82.17 82.17 78.01 𝐹𝑏∗ 3𝑖58.63 84.18 75.92 73.22 𝐹𝑎+ 3𝑖21.95 6.82 2.80 1.23 𝐹𝑏+ 3𝑖23.99 6.08 1.62 0.78 B&P 0.00 0.03 0.02 0.00 with |𝑉|≤30 within 30 minutes of CPU time. However, due to the complexity of solving the root node of the set partitioning formulation, 30 min may become not enough time to get decent results when 𝑛 increases, as we will see below. To observe a better performance of the B&P procedure over networks with 𝑛= 50 nodes, it would be necessary to increase the time limit. Moreover, we can notice that for the largest instances with 𝑛= 50 formulation (𝐹𝑏+ 3𝑖) and B&P algorithm cannot even solve the continuous relaxations within the time limit which translates into worse gap values. On the contrary, (𝐹𝑎+ 3𝑖) is able to solve its simpler continuous relaxation in most of the cases except when 𝑛= 50 and 𝑘∈ {6,8}. Finally, (𝐹𝑎∗ 3𝑖) and (𝐹𝑏∗ 3𝑖) are always able to provide lower bounds, although with poor quality, as we observe in previous sections. The results for 𝑛= 40 are really similar to the ones obtained for 𝑛= 20,30. For small values of 𝑘the best formulation is (𝐹𝑎+ 3𝑖) whereas for higher values it is the B&P algorithm. We can also remark that the performance of the B&P algorithm improves with 𝑘: one can observe that for 𝑛= 40 and 𝑘∈ {6,8}, the B&P algorithm obtains excellent results in less than 30 minutes. 4.3. Comparing clustering methods The main contribution of the paper is the introduction of a wide variety of methodologies to solve exactly the MkWNCP. Because of that, this subsection compares our exact solutions and results with the ones provided by other related (but no the same) community detection
European Journal of Operational Research 316 (2024) 519–538 532 D. Ponce et al. Table 5 Computational results of the best performance exact methods for medium-sized and largest-sized instances. 𝑛 20 30 40 50 Formulation 𝑘 24682468246 82 4 68 𝐹𝑎∗ 3𝑖Time (s) 77.00 1800.00 1800.00 1800.00 1800.00 1800.00 1800.00 1800.00 1800.00 1800.00 1800.00 1800.00 1800.00 1800.00 1800.00 1800.00 GAP (%) 0.00 39.80 56.40 49.30 38.90 66.80 70.40 76.20 60.30 79.60 82.20 82.50 69.80 85.60 87.40 87.10 Nodes 114 880 965 097 33 120 10 799 56 582 145 360 115 532 2682 215 410 29 703 10 730 2300 96 979 10 572 3031 511 Opt (%) 100 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 Dual Bound (%) 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 𝐹𝑏∗ 3𝑖Time (s) 182.00 1800.00 1800.00 1800.00 1800.00 1800.00 1800.00 1800.00 1800.00 1800.00 1800.00 1800.00 1800.00 1800.00 1800.00 1800.00 GAP (%) 0.00 47.70 34.10 19.30 50.50 68.80 76.50 66.60 69.50 83.90 85.00 82.70 81.70 89.00 90.00 89.90 Nodes 101 340 41 563 17 848 32 072 86 550 83 907 467 386 41 283 4275 452 209 7551 505 88 21 Opt (%) 100 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 Dual Bound (%) 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 𝐹𝑎+ 3𝑖Time (s) 1.90 21.42 18.22 10.50 121.97 1800.00 1800.00 1787.00 1800.00 1800.00 1800.00 1800.00 1625.00 1780.00 1800.00 1800.00 GAP (%) 0.00 0.00 0.00 0.00 0.00 4.10 1.90 0.70 16.60 11.30 5.90 3.50 23.00 24.10 70.63 92.81 Nodes 361 1928 921 214 12 318 36 096 16 876 14 966 18 413 433 318 277 1194 3 2 1 Opt (%) 100 100 100 100 100 0 0 2 0 0 0 0 20 2 0 0 Dual Bound (%) 100 100 100 100 100 100 100 100 100 100 100 100 100 96 34 8 𝐹𝑏+ 3𝑖Time (s) 22.28 45.25 11.63 2.53 1800.00 1800.00 1800.00 1177.00 1800.00 1800.00 1800.00 1800.00 1800.00 1800.00 1800.00 1800.00 GAP (%) 0.00 0.00 0.00 0.00 16.00 5.50 1.50 0.10 100.00 98.20 25.70 2.00 100.00 100.00 100.00 100.00 Nodes 304 892 98 8 1618 924 1639 3384 1 1 9 34 1 1 1 1 Opt (%) 100 100 100 100 0 0 0 66 0 0 0 0 0 0 0 0 Dual Bound (%) 100 100 100 100 100 100 100 100 0 2 78 100 0 0 0 0 B&P Time (s) 11.19 5.36 3.40 1.87 384.00 223.00 63.70 26.57 1800.00 1800.00 1222.00 608.00 1800.00 1800.00 1800.00 1800.00 GAP (%) 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 100.00 81.60 0.04 0.00 100.00 100.00 100.00 100.00 Nodes 1 2 3 3 2 7 9 6 1 2 13 23 1 1 1 1 Opt (%) 100 100 100 100 100 100 100 100 0 0 68 94 0 0 0 0 Dual Bound (%) 100 100 100 100 100 100 100 100 0 20 100 100 0 0 0 0 Table 6 Comparison between heuristic and exact solutions. 𝑛 𝑘 2 3 4 5 6 7 8 9 10 GAP_heur (%) 1.65 1.77 1.73 1.60 1.20 1.02 0.73 0.47 Opt (%) 56 47 35 38 52 59 73 90 15 GAP_heur (%) 1.78 1.38 1.53 1.67 1.73 1.68 1.58 1.33 Opt (%) 51 42 39 23 27 17 20 32 20 GAP_heur (%) 1.37 1.73 1.95 1.87 Opt (%) 44 26 16 6 30 GAP_heur (%) 1.09 1.62 1.84 1.79 Opt (%) 56 20 6 4 methods in the literature, in order to show the advantages of being able to obtain the exact optimal solution of MkWNCP. In particular, we will show the differences between the exact optimum of the MkWNCP and its approximation obtained by the heuristic spectral algorithm presented in Shi and Malik (2000), and between the exact solutions of the MkWNCP and other well-known clustering optimization models as the Maximum Modularity Problem for 𝑘-partitions. Concerning the latter, there are multiple valid MILP formulations to solve exactly. See, for instance, Agarwal and Kempe (2008), Ales and Knippel (2020), Baghersad et al. (2023), Benati and Puerto (2023), Benati et al. (2017). As we are simply comparing the solution quality and not the performances or solution times of the different methods, it is not important or necessary to point out which formulation has been used to solve exactly the MkWNCP and the Maximum Modularity Problem. Firstly, it is interesting to compare the approximation provided by the heuristic algorithm in Shi and Malik (2000) and the exact optimum over the instances that have been exactly solved (𝑛∈ {10,15,20,30}). The GAP between solutions (100(heuristic solution−optimal value)/ heuristic solution) as GAP_heur (%) and the percentage of the instances exactly solved by the heuristic (Opt (%)) are reported in Table 6. One can observe that the difference between the objective values of each of the methods is not significant because of the small GAP, less than 2%. However, as the problem becomes more complex the percentage of instances exactly solved by the heuristic is heavily reduced: for 𝑛≥20 and 𝑘≥6this percentage is less than 20%. So, we can conclude that the methodology proposed in this paper is really useful and necessary to find the exact optimum of the MkWNCP since it is not possible to do it with the current methods and algorithms that can be found in literature. Now that we are able to obtain the exact optimum of the MkWNCP thanks to our methodology, it would be interesting to test the goodness of the Normalized Cut as a community structure quality measure comparing it with other well-known objective functions in literature. In particular, we choose the Modularity to be compared with the Normalized Cut as community detection methods. As we try to compare solutions from different optimization problems, it makes little sense to compare objective values that do not represent the same function. Due to this, we will use the Normalized Mutual Information index (NMI) as a measure of 𝑘-partitions comparison, which has been applied in many works as Benati et al. (2023), Danon et al. (2005), Lancichinetti et al. (2008). The NMI metric ranges between 0 and 1, with values closer to 1 indicating strong correspondence between the 𝑘-partitions. The average NMI for each combination of parameters is reported in Table 7. In general, both partitions are really similar because of the high NMI values. However, we can observe that as the problems become harder, for intermediate 𝑘and large 𝑛, the solutions start to differ significantly. We can conclude that the Normalized Cut function is an adequate index to detect communities due to its similar performance as compared with the modularity function, one of the most studied and analyzed quality measure for community structures. The reader may note that in Tables 6 and 7for values of 𝑛= 20,30, we only report the results for 𝑘= 2,4,6,8since in our experiments we do not consider other parameter values as explained above. 4.4. Image segmentation In this subsection, we segment some real images that appear in Shi and Malik (2000) using our exact models for the 𝑘-way normalized cut.
European Journal of Operational Research 316 (2024) 519–538 533 D. Ponce et al. Table 7 NMI between exact optima of Modularity and Normalized Cut Problems. 𝑛 𝑘 23456789 10 0.85 0.80 0.75 0.77 0.83 0.88 0.92 0.99 15 0.81 0.81 0.74 0.69 0.73 0.74 0.77 0.81 20 0.78 0.68 0.64 0.69 30 0.54 0.62 0.55 0.57 The way to split these images is to consider each pixel of an image as a node of the graph and define a family of weights that measure the similarity between nodes based on the characteristics of the pixels. All the images considered in this subsection are grayscale images, so we will use the following similarity weights for any 𝑖≠𝑗: 𝑤𝑖𝑗 =𝑒 −‖𝐹(𝑖)−𝐹(𝑗)‖2 2 𝜎𝐼⋅⎧ ⎪ ⎨ ⎪ ⎩ 𝑒 −‖𝑋(𝑖)−𝑋(𝑗)‖2 2 𝜎𝑋,if ‖𝑋(𝑖) − 𝑋(𝑗)‖2< 𝑟, 0,otherwise, (79) proposed by Shi and Malik (2000), where 𝑋(𝑖)is the spatial location of node 𝑖,𝐹(𝑖)is the brightness intensity of pixel 𝑖in the grayscale image, and 𝜎𝑋,𝜎𝐼and 𝑟are a given parameter settings. 𝜎𝑋and 𝜎𝐼 can be understood as the standard deviations of the spatial location and brightness intensity of the pixels, respectively. Note that the weight 𝑤𝑖𝑗 = 0 for any pair of nodes 𝑖and 𝑗that are more than 𝑟pixels apart. Clearly, the similarity between two pixels depends on the proximity of their spatial position and their brightness intensity values. As we have noticed in the above subsection, the proposed formulations are not able to obtain the exact optimum of larger instances, so we can apply two different preprocessing techniques to the images, as a heuristic procedure, in order to reduce the problem dimension and be able to apply the different exact methods presented in the previous section. Firstly, the number of pixels can be resized and reduce the number of nodes but worsening the quality of the image. And, secondly, some kind of pre-clustering method can be used to fix the values of some variables of the formulation or to impose subset of nodes that must be in the same community. In particular, we have used the wellknown clustering method Simple Linear Iterative Clustering (SLIC) based on 𝑘-means that was first proposed by Achanta et al. (2012). This algorithm groups small subsets of pixels into superpixels reducing the dimension of the image (Ng et al.,2023). The concept of superpixel is explained and presented by Ren and Malik (2003). We have used the Python package sckit-image found in https://scikitimage.org to apply the SLIC algorithm to the images we want to segment. Then, it will be considered that each pair of pixels belonging to the same superpixel are in the same community of the final 𝑘partition. We call this new problem as the minimum 𝑘-way normalized cut under the SLIC algorithm conditions. Obviously, this new problem is not equivalent to the (MkWNCP), but it can be used as a heuristic method. The SLIC algorithm conditions can be imposed on our exact models in two ways: •Fixing the value of the variables that relate pixels to each other. •Considering each superpixel as a node of a new weighted graph and for each pair of superpixels generating a weight equal to the sum of all the weights between the pair of superpixels. In this case, loops with non-null weights can appear even if the original graph has no loops. Both processes are equivalent but the one that reduces the most the complexity of the problem is the second one. The weights between superpixels of the graph obtained by this process are called 𝑤𝑠 𝑖𝑗 . When this new weighted graph is built, to solve exactly the minimum 𝑘way normalized cut under the SLIC algorithm conditions, we must consider the weights between pixels of the same superpixel. So, given Table 8 Normalized cut value comparison between heuristic of Shi and Malik (2000) and our heuristic. Dataset Algorithm Our heuristic Shi and Malik (2000) heuristic Fig. 5 0.0192949 0.0206377 Fig. 6 0.0276930 0.0363716 Fig. 7 0.0532918 0.0595078 Fig. 8 0.0702242 0.1138075 Fig. 9 0.1876952 0.1736982 Fig. B.10a0.2150427 0.2235074 Fig. B.11b0.1732828 0.1965715 aSee Appendix B. bSee Appendix B. two superpixels 𝑖and 𝑗, that can be the same, the weight associated to the edge (𝑖, 𝑗)is equal to: 𝑤𝑠 𝑖𝑗 =∑ 𝑖′∈𝑠𝑖,𝑗′∈𝑠𝑗 𝑤𝑖′𝑗′, where 𝑠𝑖is the subset of pixels that belong to the superpixel 𝑖. This preprocess allows us to obtain a lower-sized graph whose minimum 𝑘-way normalized cut is equal to the minimum 𝑘-way normalized cut of the original graph under the conditions imposed by the SLIC preclustering. It is important to remark that we impose that the initial SLIC pre-clustering never constructs a number of superpixels strictly lower than 𝑘so that the minimum 𝑘-way normalized cut under the SLIC algorithm conditions is feasible. Finally, we have applied this heuristic process to seven images that appear in Shi and Malik (2000) using the same parameters used in that paper. For each image, we show the original image, the image obtained after applying the SLIC algorithm to obtain approximately 50 superpixels and the 𝑘different subimages into which the original image is segmented. The image segmentations are shown in Figs. 5–9and Figs. B.10– B.11. Although in Shi and Malik (2000) some of these images are treated as images of different types, we have considered that all of them are grayscale images that can be also called brightness images. We conclude that in most of the cases the subimages obtained by our heuristic process are really similar to the ones reported in Shi and Malik (2000) and the quality of our segmentation is at least equivalent or better than the one obtained by the method in Shi and Malik (2000). These results validate our methodology as an interesting approach to use normalized cuts for image segmentation. As in previous experimental subsections, the heuristic algorithm proposed in Shi and Malik (2000) and our heuristic for image segmentation have been compared to show their performances. The normalized cut objective values of each algorithm are reported in Table 8. In most of images that we have tested, the heuristic algorithm proposed in this paper outperforms the original heuristic presented in Shi and Malik (2000). 5. Conclusions In this paper, several solution methods have been introduced to solve the minimum 𝑘-way normalized cut problem. Although the benefits of solving this problem in community detection (rather than other methods, as the minimum cut problem) were known, the lack of exact algorithms to solve it, has made this technique less used than its alternatives. Here, we propose for the first time a mixed-integer linear reformulation of the original non-linear normalized cut function and analyze theoretically and empirically some MILP formulations. Furthermore, another exact method has been developed, namely a set partitioning formulation. We assess the usefulness of all these methods on an extensive computational experience both with randomly
European Journal of Operational Research 316 (2024) 519–538 534 D. Ponce et al. Fig. 5. 102 ×104 Synthetic image segmentation with parameters 𝑘= 2,𝜎𝐼= 0.2,𝜎𝑋= 5 and 𝑟= 3. Fig. 6. 106 ×106 Synthetic image segmentation with parameters 𝑘= 3,𝜎𝐼= 0.1,𝜎𝑋= 5 and 𝑟= 3. Fig. 7. 108 ×78 Planes image segmentation with parameters 𝑘= 4,𝜎𝐼= 0.01,𝜎𝑋= 4 and 𝑟= 5. Fig. 8. 122 ×89 Person image segmentation with parameters 𝑘= 5,𝜎𝐼= 0.1,𝜎𝑋= 4 and 𝑟= 5. generated instances and with a set of actual images taken from the specialized literature of image segmentation. The results reported in the paper validate our proposal in terms of solution quality and running times. Further research on this topic includes exploring alternative matheuristics and metaheuristics that help in solving faster and larger instances of these problems. Even with the limitations of complexity in the use of exact models, a large number of applications can be found in several disciplines where solution time is not essential and quality is. For a deep analysis in many areas as biology, ecology, economics, or sociology, it is more crucial to obtain a good solution in a reasonable time than an instantaneous bad solution. A clarifying example can be seen in Calvino et al. (2022), where each scanning-transmission electron microscopy (STEM) image takes several hours to be processed and the main concern is not the processing time but to get a high quality solution. It is in those cases that algorithms such as the ones proposed in this paper become really important.
European Journal of Operational Research 316 (2024) 519–538 535 D. Ponce et al. Fig. 9. 108 ×121 Weather radar image segmentation with parameters 𝑘= 6,𝜎𝐼= 0.007,𝜎𝑋= 15 and 𝑟= 10. Acknowledgments The authors of this research acknowledge financial support by: the research project PID2020-114594GB-C21 funded by the Spanish Ministerio de Ciencia e Innovación and Agencia Estatal de Investigación, Spain (MCIN/ AEI /10.13039/501100011033); and Junta de Andalucía P18-FR-1422. The first and second authors also acknowledge partial support from: NetmeetData: Ayudas Fundación BBVA a equipos de investigación científica 2019; B-FQM-322-UGR20; AT21_00032; and FQM-331. The first author also acknowledges: PAIDI 2020, postdoctoral fellowship financed by European Social Fund and Junta de Andalucía, Spain; and project CIGE/2021/161 from Generalitat Valenciana, Spain. The third author also acknowledges PREDOC PAIDI2020, predoctoral fellowship financed by European Social Fund and Junta de Andalucía, Spain. Appendix A. Branch-and-price details Here, we explain in detail some aspects of the implementation of the branch-and-price for a better understanding of the reader. A.1. Farkas pricing By Duality Theorems and Farkas’ Lemma, if 𝐹𝐷𝑃 is unbounded then 𝐹𝑅𝑆𝑃 is infeasible. In column generation procedures, this infeasibility could be avoided adding variables to the pool of columns. Hence, when 𝐹𝑅𝑆𝑃 is infeasible, the program returns Farkas values (𝛾∗, 𝜂∗)associated with constraints (66) and (67). Since in this stage we are solving a feasibility problem instead of an optimality one, primal an dual can be simplified considering 𝑐𝑆= 0. Thus, {𝑦𝑆∶ − ∑𝑖∈𝑆𝛾∗ 𝑖−𝜂∗<0} are the candidates to be added to the restricted problem in order to recover feasibility and 𝑆is built using a simplified version of pricing problem (𝐹𝑃 𝑟𝑖𝑐𝑖𝑛𝑔 ): (𝐹𝐹 𝑎𝑟𝑘𝑎𝑠𝑃 𝑟𝑖𝑐𝑖𝑛𝑔 )min − 𝑛 ∑ 𝑖=1 𝛾∗ 𝑖𝑥𝑖−𝜂∗(A.1) 𝑠.𝑡. ∶𝑥𝑖∈ {0,1},∀𝑖= 1,…, 𝑛. (A.2) Note that this problem, at the root node, can be solved easily by inspection. Additionally, when the value of the objective function is non-negative, the current node is cut off by infeasibility. Otherwise, 𝑆= {𝑖∶𝑥𝑖= 1, 𝑖 = 1,…, 𝑛}is added to Sthat is updated and 𝐹𝑅𝑆𝑃 has to be solved over. A.2. Branching When 𝐹𝑅𝑆𝑃 is solved, but the solution is not integer, the next step is to define an adequate branching rule to explore the searching tree. In this problem, we apply the Ryan and Foster branching rule (Ryan & Foster,1981) which is known to be balanced for set partitioning formulations. Given a solution with fractional 𝑦-variables in a node, it might occur that 0<∑ 𝑆∈S∶𝑖1,𝑖2∈𝑆 𝑦𝑆<1, for some 𝑖1, 𝑖2∈𝐼, 𝑖1< 𝑖2.(A.3) The Ryan and Foster branching is based in the following result: Proposition 12. If ∑𝑆∈S 𝑖1,𝑖2∈𝑆 𝑦𝑆∈ {0,1}, for all 𝑖1, 𝑖2∈𝐼, such that 𝑖1< 𝑖2, then 𝑦𝑆∈ {0,1} for all 𝑆∈S. The reader can find the proof in Barnhart et al. (1998). Provided that (A.3) happens, in order to find an integer solution, we create the following branches from the current node: left branch, where 𝑖1and 𝑖2must belong to different communities; and right branch, where 𝑖1and 𝑖2must belong to the same community. Since, in the left branch, ∑𝑆∈S 𝑖1,𝑖2∈𝑆 𝑦𝑆= 0, we have to set locally to zero all variables 𝑦𝑆 such that 𝑖1, 𝑖2∈𝑆. We have to do it for the pool of variables we have at the moment that the node is created, but also when new columns are added. Similarly, in the right branch ∑𝑆∈S 𝑖1,𝑖2∈𝑆 𝑦𝑆= 1. Thus, locally, we set to zero 𝑦𝑆such that 𝑖1∈𝑆and 𝑖2∉𝑆, or 𝑖1∉𝑆and 𝑖2∈𝑆. The reader may note that tree exploration strategy is directly done by default using SCIP implementation. The above information is easily translated to the (Farkas) pricing problem adding one constraint to each one of the branches: (1) 𝑥𝑖1+ 𝑥𝑖2≤1for the left branch; and (2) 𝑥𝑖1=𝑥𝑖2for the right branch. In practice, we could have several pairs to be candidates for branching, so we need to define a rule to decide the pair to branch. Let 𝜍𝑖1𝑖2 be the value of the sum of 𝑦variables containing nodes 𝑖1and 𝑖2in the community, i.e., 𝜍𝑖1𝑖2=∑ 𝑆∈S 𝑖1,𝑖2∈𝑆 𝑦𝑆. For instance, 𝑎𝑟𝑔 max(𝑖1,𝑖2)min{1 − 𝜍𝑖1𝑖2, 𝜍𝑖1𝑖2}, would provide pair (𝑖1, 𝑖2)to define zeroand one-branches following the most fractional branching rule. In order to test different branching strategies, we have considered multiple ways to apply the Ryan and Foster branching.
European Journal of Operational Research 316 (2024) 519–538 536 D. Ponce et al. Table A.9 Ryan and Foster test for medium-sized instances. 𝑛 20 30 𝜂 𝑘 2 4 6 8 2 4 6 8 0.5 Time (s) 11.19 5.36 3.40 1.87 384.00 223.00 63.70 26.57 Nodes 1 2 3 3 2 7 9 6 0.4 Time (s) 11.04 5.48 3.49 1.83 423.47 234.09 78.04 30.63 Nodes 1 2 3 2 2 9 15 8 0.3 Time (s) 11.21 5.40 3.44 1.88 530.33 249.59 69.87 31.13 Nodes 1 2 3 3 3 9 12 9 0.2 Time (s) 10.36 5.51 3.45 1.88 359.23 269.41 78.20 34.13 Nodes 1 2 3 3 2 10 14 10 0.1 Time (s) 10.56 5.68 3.44 1.86 430.83 309.78 85.37 34.76 Nodes 1 2 3 2 2 12 16 11 To generalize the above idea, let define 𝑎𝑟𝑔 min(𝑖1,𝑖2){|𝜂− min{1 − 𝜍𝑖1𝑖2, 𝜍𝑖1𝑖2}|}being 𝜂∈ [0,0.5]. The bigger is 𝜂 the more similar the branching selection pair is to the most fractional rule, whereas the smaller it is the selected pair furthers the selection of pairs that are about to be in same/different communities. Instead of branching on the pair with value closest to 0.5, it is possible to change the parameter for which we will branch. In particular, we have solved the instances with size 𝑛∈ {20,30} from Section 4for different values of this branching parameter (𝜂∈ {0.1,0.2,0.3,0.4,0.5}). These results are reported in Table A.9. We can conclude that the branching that selects the most fractional pair outperforms the other pair selection rules tested. Hence, we have used this rule, i.e., 𝜂= 0.5for the experiments in Section 4. A.3. Branch-and-price procedure The B&P process is represented in Algorithm 1. For the sake of readability, we detail the following functions: •Solve() has two arguments: the set Swhich define the variables 𝑦𝑆, 𝑆 ∈S, to solve 𝐹𝑅𝑆𝑃 ; and node which let the program set the local bounds according with the branching information. This function solves 𝐹𝑅𝑆𝑃 and returns feasible being (𝛾∗, 𝜂∗, 𝑦∗)the dual multipliers and the primal solution in case that the restricted relaxed master problem is feasible. If not, the function returns infeasible being (𝛾∗, 𝜂∗)the Farkas values. •SolvePricingProblem() (SolveFarkasPricingProblem()) solves 𝐹𝑃 𝑟𝑖𝑐𝑖𝑛𝑔 (𝐹𝐹 𝑎𝑟𝑘𝑎𝑠𝑃 𝑟𝑖𝑐𝑖𝑛𝑔 ) with dual multipliers (Farkas values) (𝛾∗, 𝜂∗) as data. During the solution process, more than one feasible solution could be returned. Those solutions can be sorted and create a new set 𝑆to be included in S′as long as its value of the objective function of the subproblem is negative. In this context, this is called multiple pricing. •UpdateBound() returns true if lower and upper bound coincide. We distinguish three cases represented in the third argument: 1. If the last 𝐹𝑅𝑆𝑃 solved is feasible and 𝑦∗is integer, the node is pruned by optimality and the upper bound is updated. 2. If the last 𝐹𝑅𝑆𝑃 solved is feasible and 𝑦∗is fractional, the lower bound is updated. In case this node is not pruned by bound, branching is necessary. 3. If the last 𝐹𝑅𝑆𝑃 solved is infeasible, the node is pruned by infeasibility (cutoff ). •Branch() branches fractional solution 𝑦∗according with the Ryan and Foster branching described above. Regarding variable selection, among the fractional candidates, we select the variable with value closest to 0.5, i.e., the most infeasible branching is used. Algorithm 1: pseudocode of the B&P for MkWNCP Input: An instance of MkWNCP with data 𝐺= (𝑉 , 𝐸, 𝑤). Output: An optimal partition 𝛱of 𝑉. 1S= ∅ 2optimality ←false 3node ←root node 4while optimality = false do 5(𝛾∗, 𝜂∗, 𝑦∗)←Solve(S, node) 6if Solve(S, node)is feasible then 7S′←SolvePricingProblem(𝛾∗, 𝜂∗, node) 8else 9S′←SolveFarkasPricingProblem(𝛾∗, 𝜂∗, node) 10 if S′≠∅then 11 S←S∪S′ 12 else 13 if 𝑦∗is integer then 14 if UpdateBound(𝑦∗, node, integer) then 15 optimality ←true 16 else 17 node ←NextNode(𝑀𝑃 ) 18 else if 𝑦∗is fractional then 19 if UpdateBound(𝑦∗, node, fractional) then 20 optimality ←true 21 else 22 if node is not pruned then 23 Branch(𝑦∗) 24 node ←NextNode(𝑀𝑃 ) 25 else 26 if UpdateBound(NULL, node, cutoff) then 27 optimality ←true 28 else 29 node ←NextNode(𝑀𝑃 ) •Nextnode() returns the next node to be processed among the candidates in order to solve the master problem. Appendix B. Extra image segmentations See Fig. B.10 and B.11.
European Journal of Operational Research 316 (2024) 519–538 537 D. Ponce et al. Fig. B.10. 115 ×79 Baseball game image segmentation with parameters 𝑘= 7,𝜎𝐼= 0.1,𝜎𝑋= 4 and 𝑟= 5. Fig. B.11. 132 ×81 Zebra image segmentation with parameters 𝑘= 7,𝜎𝐼= 0.1,𝜎𝑋= 4 and 𝑟= 5. References Achanta, R., Shaji, A., Smith, K., Lucchi, A., Fua, P., & Süsstrunk, S. (2012). Slic superpixels compared to state-of-the-art superpixel methods. IEEE Transactions on Pattern Analysis and Machine Intelligence,34(11), 2274–2282. Agarwal, G., & Kempe, D. (2008). Modularity-maximizing graph communities via mathematical programming. The European Physical Journal B,66, 409–418. Ales, Z., & Knippel, A. (2020). The k-partitioning problem: Formulations and branch-and-cut. Networks,76(3), 323–349. Aloise, D., Cafieri, S., Caporossi, G., Hansen, P., Perron, S., & Liberti, L. (2010). Column generation algorithms for exact modularity maximization in networks. Physical Review. E, Statistical, Nonlinear, and Soft Matter Physics,82, Article 046112. Arredondo, V., Martínez-Panero, M., Peña, T., & Ricca, F. (2021). Mathematical political districting taking care of minority groups. Annals of Operations Research,305(1), 375–402. Baghersad, M., Emadikhiav, M., Huang, C. D., & Behara, R. S. (2023). Modularity maximization to design contiguous policy zones for pandemic response. European Journal of Operational Research,304(1), 99–112. Barnhart, C., Johnson, E., Nemhauser, G., Savelsbergh, M., & Vance, P. (1998). Branch-and-price: Column generation for solving huge integer programs. Operations Research,46, 316–329. Barrientos, M., & Madrid, H. (2011). Normalized cut based edge detection. In J. F. Martínez-Trinidad, J. A. Carrasco-Ochoa, C. Ben-Youssef Brants, & E. R. Hancock (Eds.), Pattern recognition (pp. 211–219). Berlin, Heidelberg: Springer Berlin Heidelberg, ISBN: 978-3-642-21587-2. Benati, S., Ponce, D., Puerto, J., & Rodríguez-Chía, A. M. (2022). A branch-andprice procedure for clustering data that are graph connected. European Journal of Operational Research,297(3), 817–830. Benati, S., & Puerto, J. (2023). A network model for multiple selection questions in opinion surveys. Quality & Quantity, 1–17. Benati, S., Puerto, J., & Rodríguez-Chía, A. M. (2017). Clustering data that are graph connected. European Journal of Operational Research,261(1), 43–53. Benati, S., Puerto, J., Rodríguez-Chía, A. M., & Temprano, F. (2022). A mathematical programming approach to overlapping community detection. Physica A: Statistical Mechanics and its Applications,602, Article 127628. Benati, S., Puerto, J., Rodríguez-Chía, A. M., & Temprano, F. (2023). Overlapping communities detection through weighted graph community games. PLoS One,18(4), 1–35. http://dx.doi.org/10.1371/journal.pone.0283857. Bestuzheva, K., Besançon, M., Chen, W.-K., Chmiela, A., Donkiewicz, T., van Doornmalen, J., Eifler, L., Gaul, O., Gamrath, G., Gleixner, A., Gottwald, L., Graczyk, C., Halbig, K., Hoen, A., Hojny, C., van der Hulst, R., Koch, T., Lübbecke, M., Maher, S. J., .... Witzig, J. (2021). The SCIP optimization suite 8.0:ZIB-Report 21-41, Zuse Institute Berlin, URL http://nbn-resolving.de/urn:nbn:de:0297-zib-85309. Blanco, V., Gázquez, R., Ponce, D., & Puerto, J. (2023). A branch-and-price approach for the continuous multifacility monotone ordered median problem. European Journal of Operational Research,305(1), 105–126. Blanco, V., Japón, A., Ponce, D., & Puerto, J. (2021). On the multisource hyperplanes location problem to fitting set of points. Computers & Operations Research,128, Article 105124. Cai, W., Wu, J., & Chung, A. C. (2006). Shape-based image segmentation using normalized cuts. In 2006 international conference on image processing (pp. 1101–1104). Calvino, J. J., López-Haro, M., Muñoz-Ocaña, J. M., Puerto, J., & Rodríguez-Chía, A. M. (2022). Segmentation of scanning-transmission electron microscopy images using the ordered median problem. European Journal of Operational Research, [ISSN: 0377-2217] 302(2), 671–687. http://dx.doi.org/10.1016/j.ejor.2022.01.022. Carrizosa, E., Marín, A., & Pelegrín, M. (2020). Spotting key members in networks: Clustering-embedded eigenvector centrality. IEEE Systems Journal,14(3), 3916–3925. Cheeger, J. (1971). A lower bound for the smallest eigenvalue of the Laplacian (pp. 195–200). Princeton: Princeton University Press, ISBN: 978-1-4008-6931-2. Costa, A. (2015). MILP formulations for the modularity density maximization problem. European Journal of Operational Research,245(1), 14–21. Costa, A., Ng, T. S., & Foo, L. X. (2017). Complete mixed integer linear programming formulations for modularity density based clustering. Discrete Optimization, (25), 141–158. Danon, L., Díaz-Guilera, A., Duch, J., & Arenas, A. (2005). Comparing community structure identification. Journal of Statistical Mechanics: Theory and Experiment, 2005(09), P09008. http://dx.doi.org/10.1088/1742-5468/2005/09/P09008. Deleplanque, S., Labbé, M., Ponce, D., & Puerto, J. (2020). A branch-price-and-cut procedure for the discrete ordered median problem. INFORMS Journal on Computing, 32(3), 582–599. Desaulniers, G., Desrosiers, J., & Solomon, M. M. (2005). column generation. New York: Springer. Donath, W. E., & Hoffman, A. J. (1973). Lower bounds for the partitioning of graphs. IBM Journal of Research and Development,17(5), 420–425. Fiedler, M. (1975). A property of eigenvectors of nonnegative symmetric matrices and its applications to graph theory. Czechoslovak Mathematical Journal,25(4), 619–633. Goldschmidt, O., & Hochbaum, D. S. (1994). A polynomial algorithm for the k-cut problem for fixed k. Mathematics of Operations Research,19(1), 24–37. Hansen, P., & Jaumard, B. (1997). Cluster analysis and mathematical programming. Mathematical Programming,79(1), 191–215.
European Journal of Operational Research 316 (2024) 519–538 538 D. Ponce et al. Hansen, P., Ruiz, M., & Aloise, D. (2012). A VNS heuristic for escaping local extrema entrapment in normalized cut clustering. Pattern Recognition,45(12), 4337–4345. He, L., & Zhang, H. (2016). Iterative ensemble normalized cuts. Pattern Recognition,52, 274–286. Korf, R. E. (1998). A complete anytime algorithm for number partitioning. Artificial Intelligence, [ISSN: 0004-3702] 106(2), 181–203. http://dx.doi.org/10.1016/S00043702(98)00086-1. Lancichinetti, A., Fortunato, S., & Kertész, J. (2008). Detecting the overlapping and hierarchical community structure in complex networks. New Journal of Physics, 11(3). Lübbecke, M., & Desrosiers, J. (2005). Selected topics in column generation. Operations Research,53(6), 1007–1023. Maravalle, M., Simeone, B., & Naldini, R. (1997). Clustering on trees. Computational Statistics & Data Analysis,24(2), 217–234. Matić, D., & Grbić, M. (2020). Partitioning weighted metabolic networks into maximally balanced connected partitions. In 2020 19th international symposium INFOTEH-JAHORINA (pp. 1–6). Mehrotra, A., & Trick, M. A. (1998). Cliques and clustering: A combinatorial approach. Operations Research Letters,22(1), 1–12. Ng, T. C., Choy, S. K., Lam, S. Y., & Yu, K. W. (2023). Fuzzy superpixel-based image segmentation. Pattern Recognition,134, Article 109045. Ojeda-Ruiz, I., & Lee, Y. J. (2020). A fast constrained image segmentation algorithm. Results in Applied Mathematics,8, Article 100103. Pothen, A., Simon, H. D., & Liou, K. P. (1990). Partitioning sparse matrices with eigenvectors of graphs. SIAM Journal on Matrix Analysis and Applications,11(3), 430–452. Ren, X., & Malik, J. (2003). Learning a classification model for segmentation. In Proceedings ninth IEEE international conference on computer vision:vol. 1, (pp. 10–17). Ricca, F., Scozzari, A., & Simeone, B. (2013). Political districting: From classical models to recent approaches. Annals of Operations Research,204(1), 271–299. Ryan, D., & Foster, B. (1981). An integer programming approach to scheduling. Computer Scheduling of Public Transport, (1), 269–280. Shi, J., & Malik, J. (2000). Normalized cuts and image segmentation. IEEE Transactions on Pattern Analysis and Machine Intelligence,22(8), 888–905. Sukeda, I., Miyauchi, A., & Takeda, A. (2023). A study on modularity density maximization: Column generation acceleration and computational complexity analysis. European Journal of Operational Research, [ISSN: 0377-2217] 309(2), 516–528. http://dx.doi.org/10.1016/j.ejor.2023.01.061. Tepper, M., Musé, P., Almansa, A., & Mejail, M. (2011). Automatically finding clusters in normalized cuts. Pattern Recognition,44(7), 1372–1386. Validi, H., Buchanan, A., & Lykhovyd, E. (2022). Imposing contiguity constraints in political districting models. Operations Research,70(2), 867–892. Vincken, S. (2021). computation of normalized cuts of graphs (Master thesis), Utrecht University. Wang, C., Chen, X., Nie, F., & Huang, J. Z. (2022). Directly solving normalized cut for multi-view data. Pattern Recognition,130, Article 108809. Wertheimer, M. (1938). Laws of organization perceptual forms. In W. Ellis (Ed.), Source book of gestalt psychology. Brace: Harcourt. Xu, L., Li, W., & Schuurmans, D. (2009). Fast normalized cut with linear constraints. In IEEE conference on computer vision and pattern recognition (pp. 2866–2873). Yang, J., Yang, X., Zhou, Z.-B., & Liu, Z.-Y. (2022). Graph matching based on fast normalized cut and multiplicative update mapping. Pattern Recognition,122, Article 108228. Zhong, G., & Pun, C. M. (2022). Improved normalized cut for multi-view clustering. IEEE Transactions on Pattern Analysis and Machine Intelligence,44(12), 10244–10251. Zhou, X., Wang, H., Ding, B., Hu, T., & Shang, S. (2019). Balanced connected task allocations for multi-robot systems: An exact flow-based integer program and an approximate tree-based genetic algorithm. Expert Systems with Applications,116, 10–20.
Chapter 9 Modularity for Hypergraph Clustering: Methodologies and Applications 121
(2020). A consistent option, followed by this paper, is ∆(e) = |e|(|e|−1) 2, since |e|(|e|−1) 2edges are generated for each hyperedge e. Finally, the modularity function (1) developed for simple graphs can be extended to hypergraphs with the same formula applied to the induced graph. Consider a node partition {C1, . . . , Cq}, then its modularity, later on defined as projected modularity, is: Q1=1 P e∈E |e|(|e|−1) ∆(e) q X s=1 X i,j∈Cs X e∈E Ae ij ∆(e)− P e∈E i∈e |e|−1 ∆(e) P e∈E j∈e |e|−1 ∆(e) P e∈E |e|(|e|−1) ∆(e) .(3) 3.2. Hypergraph modularity by inclusiveness The previous approach has the merit of its simplicity. However, the method is not flawless. The main ambiguity is the lack of the definition of internal and external arcs, a concept that is at the root of the Equation (1). In a simple graph, an arc is within a community if its nodes belong to the same cluster, otherwise it is outside of the communities. In a simple graph, this is a true or false condition, given that a pair belongs or does not belong to the same community. However, when applied to hypergraphs, this notion is vanishing, and actually, it does not appear in Equation (3). Indeed, in that formula a hyperarc can be internal to more than one community and being external too. For example, if e={v1, v2, v3, v4} and supposing that in modularity optimization v1and v2 are in one cluster, v3and v4are in a second cluster, then the hyperarc e, through the arcs (v1, v2)and (v3, v4) is internal to both clusters, but in a sense, it is also external because all simple arcs between nodes of different clusters, such as (v1, v3), are external. Clearly, if one wants to take advantage of algorithms developed for graphs to hypergraphs as well, then this ambiguity can be forgotten. However, one may also wonder whether a new definition of modularity can be defined, based on a clear dichotomy between internal and external hyperarcs. In the following, we are introducing a new measure of hypergraph modularity. We retain the main logical steps that were used to define Equation (1). They are: First, to define whether a hyperarc is internal or not to some cluster; second, to compare its presence (or absence) to its expected value, given that the hypergraph emerged as the random occurrence of the configuration model. Moreover, we extend that model by considering that a hyperarc can belong to a cluster with different level of membership, as it can be the case that only a fraction of its vertices are contained in a cluster. Therefore, if it is established that the hyperarc is internal or not to a cluster, then it must be weighted by its membership level. Step 1: Defining internal hyperedges: We define a hyperedge as internal to subset S⊆V, if it is incident to at least a given fraction of Snodes. More formally, let γ∈[0,1] be the contribution parameter, hyperedge e is internal to S, if and only if |S∩e| ≥ γ|e|. Analogously, for each hyperedge e, we can define a lower bound γe=⌈γ|e|⌉ such that eis internal to S, if and only if |S∩e| ≥ γe. Relevant values used in our contribution are γ > 0.5, but values below that threshold are still admissible. Let Ae Sbe defined as: Ae S= 1,if hyperedge eis internal of S, i.e., |S∩e| ≥ γe, 0,otherwise. The number ASof the internal hyperedges of Sis: AS=|{e∈E:|S∩e| ≥ γe}|, so that AS=P e∈E Ae S. Step 2: The weight of a hyperedge: Hyperedges should have different weights depending on the fraction of their nodes that are contained in a community. For example, consider Figure 3, there, a cluster Cs⊆Vof four nodes is coloured in orange, and a fifth node not belonging to Csis coloured in green. Two hyperedges are drawn. The first e1is incident to all the orange nodes, the second e2in incident to two orange and to the green nodes. For γ= 0.5, both e1and e2are internal to the community Cs, however, we should take into account that e1is more inside than e2, and therefore it should be more important than e2to define Cs. Therefore, weights 6
w(e, Cs)should be defined to measure the number of nodes of ethat are contained in Cs. If eis internal to Cs, we propose to use the value w(e, Cs)=(|Cs∩e|−γe+1), being γe:= ⌈γ|e|⌉, which corresponds to the number of internal nodes that exceed the threshold γe. If eis not internal to Cs, then w(e, Cs) = 0. So, these weights can be defined as w(e, Cs) = (|Cs∩e| − γe+ 1)Ae Cs. Note that if the modularity is defined on graphs, there is no need to consider weights: in that case w(e, Cs)=1or 0, depending on whether e= (i, j)⊆Csor not. e1e2 γ= 0.5 Figure 3: e1and e2are internal hyperedges of orange community for γ= 0.5. Step 3: The configuration model for hypergarphs: Given a vertex clustering, modularity is an index that compares the number of the arcs internal to clusters with their expected value conditional to the configuration model. In the following, we are developing formulas for the expected number of arcs that are internal to partitions for the case of hypergraphs. The standard configuration model is a random graph model, in which the vertex degrees are the same, but arcs are rewired in such a way that the resulting graph has no structure. The configuration model of the hypergraph Hproduces random occurrences of a hypergraph with the same node degrees of H, the same number of hyperedges as H, finally, all random hyperedges have the same size as the hyperedges of H. A simple way to understand the random generation process of the configuration model is to represent each hyperedge as a facility to which the nodes that belong to that hyperedge are connected, see Figure 4. So, the hypergraph H= (V, E)is represented by a bipartite graph B= (V, E, Q), in which one vertex class Vrepresents the nodes, the other vertex class Erepresents the edges, and there is an arc (v, e)∈Qif v∈e. 1 3 2 4 5 e1e2e3 e4 Figure 4: Alternative representation of a hypergraph from Figure 2. As in the original configuration model, every arc of Bis cut into two stubs, one stub incident to a facility (hyperarc) node, the other stub incident to a node. Therefore, we obtain δnode stubs and δfacility stubs. All the random graphs of the configuration model are all the matchings between the two families of stubs. In Figure 5, we have drawn the process of obtaining a random configuration from a given occurrence of H. 7
12 43 5 e3e1e2 e4 12 43 5 e3e1e2 e4 12 43 5 e3e1e2 e4 Figure 5: New configuration model process example with hypergraph from figures 2 and 4. Since every random hypergraph is defined by a given matching between δnode stubs and δfacility stubs, then there are δ!random hypergraphs that can be obtained by the configuration model, all with the same probability. Consider a subset S⊆V, the adjacency degree of Sis δS=P i∈S δi. Next, consider Sand a hyperarc e. In some random hypergraphs, econtains exactly tnodes of S. That is, exactly tstubs incident to eare paired with telements of S. The number of hypergraphs with this property is: δS tδ−δS |e| − t|e|!(δ− |e|)!. Indeed, there are δS tpossible subsets of tstubs that can be selected out of Sand δ−δS |e|−tpossible subsets of |e| − tstubs can be selected from V\S. Those |e|stubs can be matched in |e|!possible ways with the stubs incident to e. Moreover, these combinations must be combined with the remaining (δ− |e|)! possible matchings between the rest of the stubs. Finally, if we divide that value by the total number of possible cases, the probability of the random occurrence of a hypergraph in which |S∩e|=tis: δS tδ−δS |e|−t|e|!(δ− |e|)! δ!=δS tδ−δS |e|−t δ |e|. One can observe that the probability corresponds to the hypergeometric distribution in which we are taking |e|elements from a population of δelements. There are tof these |e|elements that belong to a group of cardinality δScorresponding to the stubs that are adjacent to S. The remaining |e|−tare from a group of cardinality δ−δS, corresponding to the stubs not adjacent to S. Note that the above probability only makes sense if max{0,|e|+ δS−δ} ≤ t≤min{δS,|e|}. So, considering we need at least γeconnections between Sand e, the exact expected value of Ae Sis: E[Ae S] = min{|e|,δS} X t=max{γe,|e|+δS−δ} δS tδ−δS |e|−t δ |e|. Let us notice that the above expression is equal to the probability of ebeing internal to S. Step 4: Modularity for hypergraphs. Modularity compares the weighted sum of the internal hyperedges to community to the expected value of the same weighted sum under the configuration model. The weighted sum of the internal hyperedges to community can be expressed as: q X s=1 X e∈E (|Cs∩e| − γe+ 1)Ae Cs. Therefore, summarizing our discussion above, we get to the following hypergraph modularity equation: 8
Qγ 2=1 m q X s=1 X e∈E (|Cs∩e| − γe+ 1)Ae Cs− min{|e|,δCs} X t=max{γe,|e|+δCs−δ} (t−γe+ 1)δCs tδ−δCs |e|−t δ |e| .(4) In particular, the expression (|Cs∩e| − γe+ 1)Ae Cstakes value 0if eis not internal to community Cs, otherwise it takes the value of the weight (|Cs∩e| − γe+ 1). Then the same expression is compared to its exact expected value. The two modularity measures (3) and (4) can be used to determine the best partition of V. In the following, we will introduce two Mixed Integer Linear Programming (MILP) formulations to that purpose. 3.2.1. MILP formulations for Maximum Modularity The projected modularity (3) is maximized through the clique partition model. The following variables define the partition: xij = 1,if nodes iand jare in the same community, 0,otherwise, for any i, j = 1, . . . , n, such that i<j. The maximum modularity/clique partition model is the following integer linear program: (F1 x) : max 1 P e∈E |e|(|e|−1) ∆(e) n−1 X i=1 n X j=i+1 X e∈E Ae ij ∆(e)− P e∈E i∈e |e|−1 ∆(e) P e∈E j∈e |e|−1 ∆(e) P e∈E |e|(|e|−1) ∆(e) xij (5a) s.t. xij +xjk −xik ≤1,∀i, j, k = 1, . . . , n, i < j < k, (5b) xij −xjk +xik ≤1,∀i, j, k = 1, . . . , n, i < j < k, (5c) −xij +xjk +xik ≤1,∀i, j, k = 1, . . . , n, i < j < k, (5d) xij ∈ {0,1}, i, j = 1, . . . , n, i < j. (5e) In this formulation, it is only necessary to impose that the binary variables xij guarantee the triangle inequalities defined by (5b)-(5d). That is, if i, j ∈Csand j, k ∈Csfor some s, then we must have i, k ∈Cs. There is no need to impose additional constraints, such as the connectivity of Cs, as it has been proven that optimal partitions are composed of connected subsets in Dinh and Thai (2015). Optimizing the hypergraph modularity (4) is more complex. The potential maximum number of communities to partition the vertex set V, with |V|=n, is n. Therefore, we assume that Vis always partitioned into n communities {C1, . . . , Cn}, some of them can be empty. The partition is defined by the variables: zis = 1,if node ibelongs to community s, 0,otherwise. If Csis a non-empty community, then sis defined by the highest index of the nodes that belongs to it, i.e., if |Cs|>0then maxi∈Csi=s. So, we can define the zis variables for any i, s = 1, . . . , n, such that i≤s. In this way, symmetric solutions of a MILP problem are removed, as can be seen in Ales and Knippel (2020), Benati et al. (2022, 2017), Ponce et al. (2024). Since the hypergraph modularity (4) depends on the number of internal hyperedges and the adjacent degree 9
of each community, the following variables are needed: dls = 1,if the adjacency degree of community sis equal to l, 0,otherwise. hes = 1,if eis an internal hyperedge of community s, 0,otherwise. Variables dls are defined for any s= 1, . . . , n, and l=δs, . . . , δ −n P i=s+1 δi. Variables hes are defined for any s= 1, . . . , n and e∈E. Finally, continuous non-negative variables h+ es are needed to model the excess of an internal hyperedge with respect to γe, and they are defined for any s= 1, . . . , n and e∈E. The maximum of equation (4) can be obtained by the following MILP problem: (F2 z,γ) : max 1 m n X s=1 X e∈E hes +h+ es − δ− n P i=s+1 δi X l=δs min{|e|,l} X t=max{γe,|e|+l−δ} (t−γe+ 1)l tδ−l |e|−t δ |e| dls (6a) s.t. n X s=i zis = 1,∀i= 1, . . . , n, (6b) zis ≤zss,∀s= 1, . . . , n, i = 1, . . . , s, (6c) δ− n P i=s+1 δi X l=δs dls =zss,∀s= 1, . . . , n, (6d) δ− n P i=s+1 δi X l=δs ldls = s X i=1 δizis,∀s= 1, . . . , n, (6e) γehes +h+ es ≤ s X i=1 i∈e zis,∀s= 1, . . . , n, e ∈E, (6f) h+ es ≤hes(|e| − γe),∀s= 1, . . . , n, e ∈E, (6g) zis ∈ {0,1},∀s= 1, . . . , n, i = 1, . . . , s, (6h) dls ∈ {0,1},∀s= 1, . . . , n, l =δs, . . . , δ − n X i=s+1 δi,(6i) hes ∈ {0,1},∀s= 1, . . . , n, e ∈E, (6j) h+ es ≥0,∀s= 1, . . . , n, e ∈E. (6k) Constraints (6b)-(6c) guarantee that each node belongs to one and only one community, that a community sis non-empty if and only if the node sbelongs to it, and that sis also the maximum index of the community Cs. Then, equalities (6d)-(6e) impose that each non-empty community is assigned one and only one adjacency degree and this has to be equal to the sum of adjacency degrees of its nodes. Finally, the family of inequalities (6f) ensures that a hyperedge eis internal of a community sif and only if the number of nodes that belongs to the hyperedge and the community is greater or equal than γe. All the variables are defined as binary by (6h)-(6j). To reduce the problem size, the h-variables can be preprocessed: given a hyperedge eand a community s, if |{i∈e:i≤s}| < γe, then hes = 0 because it is impossible for eto be an internal hyperedge of community s. For the special case that γ= 1, that is, hyperedges are internal only if all their nodes are inside a community, 10
the problem (F2 z,γ)can be simplified. All the hvariables above are replaced by the following: he= 1,if eis an internal hyperedge, 0,otherwise, to obtain the formulation: (F2 z,1) : max 1 m X e∈E he− n X s=1 δ− n P i=s+1 δi X l=1 X e∈E min{|e|,l} X t=|e| l tδ−l |e|−t δ |e| dls (7a) s.t.(6b) −(6e),(6h) −(6j)., he≤1−zis,∀e∈E, s = 1, . . . , n, max i′∈ei′> s, i ∈e, i ≤s, (7b) he≤1 + zjs −zis,∀e∈E, s = 1, . . . , n, max i′∈ei′≤s, i, j ∈e, i =j. (7c) In this special case, if there are two nodes from a hyperedge ethat belong to different communities then ecan not be internal, due to constraints (7b)-(7c). 4. Computational tests In this section, we analyze how the modularity equations (3) and (4), though their optimization models, can recover the community structure of hypergraphs. The following design is common in the computational experiments in the area of network analysis, see Lancichinetti et al. (2008), Tandon et al. (2021), Benati et al. (2023). Following a rationale similar to the one for the extension of the modularity measures for hypergraphs, we have extended the procedure introduced by Lancichinetti et al. (2008), to generate hypergraphs with a given community structure using of the following parameters: •n: the number of nodes. •m: the number of hyperedges. •nc: the maximum number of communities. •L: the maximum hyperedge size. The community structure of a hypergraph can be easy or hard to recover, depending on the following parameters: •1−µ: fraction of hyperedges that are internal. •1−ξ: fraction of nodes from an internal hyperedge that belong to the same community. The noise parameter 1−ξis related to the parameter γthat is used in model (F2 z,γ)to determine an internal arc, therefore that model has been tested with γ= 1 −ξ. Moreover, model (F1 x)has been solved for all values: ∆(e) = {1,|e|,|e| − 1,|e|(|e|−1) 2}. These values include our suggestion ∆(e) = |e|(|e|−1) 2and all the other parameters suggested by the literature to reduce hypergraphs clustering to a graph problem. In the results below, we report only the best solution from graph reduction. The quality of the community structure obtained by each method is measured comparing the outcome with the true community structure, that is known by construction. The statistics that we consider are: •Normalized Mutual Information (NMI), presented in Fred and Jain (2003). •Omega index (OI), presented in Collins and Dent (1988). 11
As (F1 x)and (F2 z,γ)are MILP models formulated to solve NP-hard problems, they could imply extensive computational times to calculate the optimum. Therefore, we cannot test the models using large size hypergraphs. However, and conversely to common practices, we did not rely on heuristic procedures, but we have calculated the exact best solutions. As a consequence, our results are a reliable comparison between measures (3) and (4), because the outcomes are not biased by the quality of the heuristics. So, tests are run on instances with sizes n={10,20,30},m={50,100}, L ={3,5}and the maximum number of communities nc= 5. Tested parameter values are µ={0,0.2,0.4}and ξ={0,0.2,0.4}. Algorithms and models have been implemented in Python, and MILP formulations have been solved with Gurobi. For every combination of parameters, we have solved 50 random instances of the hypergraph modularity problem with both models. Results, reported as averages, can be see in tables 1 and 2. In most cases, according to both statistics NMI and OI, the formulation (F2 z,γ)performs better than the best partition calculated by (F1 x), with varying values of ∆(e). Indeed, when n= 10 and considered 36 parameters combination, (F2 z,γ)outperformed (F1 x)in 23 or 22 cases (depending by the index). In only one case, the reverse result has been obtained. For n= 20,(F2 z,γ)obtained the best results in 19 out of the 36 cases. For n= 30 the result has been 17 against 8. Moreover, both formulations obtained statistics close to 1, showing that hypergraph modularity is an effective concept to determine the hidden communities. Formulation F1 xF2 z,γ n m µ ξ L NMI OI NMI OI 10 50 0 031 1 1 1 51 1 1 1 0.2 31 1 1 1 5 0.98 0.96 1 1 0.4 3 0.95 0.91 0.96 0.94 5 0.96 0.92 0.98 0.94 0.2 031 1 0.98 0.96 5 0.97 0.96 1 1 0.2 3 0.98 0.97 1 1 5 0.92 0.85 0.98 0.96 0.4 3 0.81 0.67 0.90 0.81 5 0.95 0.88 0.96 0.92 0.4 030.94 0.89 0.94 0.87 5 0.87 0.80 1 1 0.2 3 0.93 0.87 0.96 0.92 5 0.83 0.73 0.92 0.87 0.4 3 0.81 0.72 0.90 0.82 5 0.64 0.48 0.84 0.70 100 0 031 1 1 1 51 1 1 1 0.2 31 1 1 1 51 1 1 1 0.4 3 0.98 0.96 1 1 5 0.98 0.97 1 1 0.2 031 1 1 1 5 0.96 0.93 1 1 0.2 31 1 1 1 51 1 1 1 0.4 3 0.96 0.93 1 1 5 0.92 0.85 1 1 0.4 031 1 1 1 5 0.94 0.90 1 1 0.2 3 0.96 0.93 0.99 0.96 5 0.92 0.87 1 1 0.4 3 0.96 0.91 0.98 0.94 5 0.85 0.75 1 1 Formulation F1 xF2 z,γ n m µ ξ L NMI OI NMI OI 20 50 0 030.97 0.96 0.97 0.95 51 1 1 1 0.2 30.95 0.92 0.95 0.91 51 1 1 1 0.4 3 0.86 0.80 0.92 0.87 5 0.98 0.98 0.99 0.99 0.2 03 0.88 0.82 0.90 0.85 51 1 0.97 0.95 0.2 3 0.92 0.87 0.95 0.93 50.96 0.93 0.94 0.88 0.4 3 0.67 0.51 0.79 0.65 5 0.91 0.87 0.94 0.91 0.4 03 0.75 0.67 0.77 0.69 5 0.96 0.95 0.97 0.96 0.2 3 0.73 0.63 0.80 0.72 50.86 0.80 0.85 0.78 0.4 3 0.43 0.25 0.46 0.27 5 0.75 0.66 0.80 0.70 100 0 031 1 1 1 51 1 1 1 0.2 31 1 1 1 51 1 1 1 0.4 3 0.96 0.94 0.98 0.97 5 0.98 0.97 0.99 0.98 0.2 031 1 1 1 51 1 1 1 0.2 30.99 0.99 0.99 0.99 51 1 1 1 0.4 3 0.89 0.84 0.92 0.89 5 0.97 0.95 1 1 0.4 03 0.93 0.90 0.94 0.91 51 1 1 1 0.2 3 0.91 0.88 0.92 0.89 51 1 1 1 0.4 3 0.61 0.48 0.65 0.48 50.88 0.82 0.88 0.85 Table 1: Computational results with n= 10,20. 12
Formulation F1 xF2 z,γ n m µ ξ L NMI OI NMI OI 30 50 0 030.91 0.86 0.90 0.87 50.99 0.98 0.99 0.98 0.2 30.95 0.92 0.94 0.92 50.98 0.96 0.96 0.93 0.4 3 0.71 0.58 0.76 0.65 5 0.89 0.85 0.96 0.93 0.2 030.83 0.76 0.80 0.72 50.98 0.97 0.98 0.97 0.2 3 0.80 0.71 0.84 0.79 5 0.89 0.85 0.92 0.87 0.4 3 0.51 0.33 0.58 0.43 50.77 0.70 0.77 0.70 0.4 03 0.51 0.36 0.54 0.42 50.73 0.67 0.72 0.65 0.2 3 0.59 0.45 0.60 0.46 5 0.75 0.66 0.80 0.75 0.4 30.38 0.20 0.38 0.23 5 0.48 0.38 0.52 0.40 100 0 030.99 0.99 0.99 0.99 51 1 1 1 0.2 30.99 0.99 0.99 0.99 51 1 1 1 0.4 3 0.90 0.88 0.94 0.92 50.99 0.99 0.98 0.98 0.2 03 0.97 0.96 0.98 0.97 51 1 1 1 0.2 30.98 0.98 0.98 0.97 50.99 0.99 0.99 0.99 0.4 3 0.69 0.61 0.74 0.65 5 0.96 0.95 0.98 0.97 0.4 030.72 0.64 0.72 0.62 5 0.97 0.97 0.99 0.99 0.2 3 0.72 0.62 0.74 0.68 5 0.95 0.95 0.96 0.95 0.4 30.47 0.35 0.45 0.28 5 0.70 0.65 0.72 0.67 Table 2: Computational results with n= 30. Grouping the instances further, therefore reporting the same results, but considering the average marginal metrics after pairing nwith m, or L, or µor ξ, results reported in Table 3 are reported. First, it is interesting to observe that the quality of the partition improves as the number of hyperedges increases. We may guess that it happens because the ideal communities are better connected, or well separated, when the number of hyperedges is greater. Something similar happens with the parameter L, because the greater L, the easier is to cover and to connect the nodes that within a community. Conversely, and as could be expected, to increase the noise parameters µand ξ, affects negatively the quality of the partitions. We outline that parameter ξhas a more negative impact on F1 x, because in that model a noising hyperarc contributes to the definition of a community. Conversely, in the second formulation F2 z,γ that noising factor is taken into account by γ, and therefore it does not negatively affects the index. As a final summary, for each size nthe overall average values are reported in Table 4. As can be seen, we can conclude that the hypergraph modularity equation (4) outperforms the projected modularity (3). 13
Formulation F1 xF2 z,γ n m NMI OI NMI OI 10 50 0.92 0.87 0.96 0.93 100 0.97 0.94 1 0.99 20 50 0.87 0.81 0.89 0.83 100 0.95 0.93 0.96 0.94 30 50 0.76 0.68 0.77 0.70 100 0.89 0.86 0.90 0.87 Formulation F1 xF2 z,γ n L NMI OI NMI OI 10 3 0.96 0.93 0.98 0.96 5 0.93 0.88 0.98 0.97 20 3 0.86 0.80 0.88 0.83 50.96 0.94 0.96 0.94 30 3 0.76 0.68 0.77 0.70 5 0.89 0.86 0.90 0.87 Formulation F1 xF2 z,γ n µ NMI OI NMI OI 10 00.99 0.98 0.99 0.99 0.2 0.96 0.92 0.98 0.97 0.4 0.89 0.82 0.96 0.92 20 00.98 0.96 0.98 0.97 0.2 0.93 0.90 0.95 0.92 0.4 0.82 0.75 0.84 0.77 30 0 0.94 0.92 0.95 0.93 0.2 0.86 0.82 0.88 0.84 0.4 0.66 0.57 0.68 0.59 Formulation F1 xF2 z,γ n ξ NMI OI NMI OI 10 0 0.97 0.96 0.99 0.99 0.2 0.96 0.93 0.99 0.98 0.4 0.90 0.83 0.96 0.92 20 00.96 0.94 0.96 0.94 0.2 0.94 0.92 0.95 0.93 0.4 0.82 0.76 0.86 0.80 30 00.88 0.85 0.88 0.85 0.2 0.88 0.84 0.89 0.86 0.4 0.71 0.62 0.73 0.65 Table 3: Marginal values for each one of the 4-factors in the experiment Formulation F1 xF2 z,γ n NMI OI NMI OI 10 0.94 0.91 0.98 0.96 20 0.91 0.87 0.92 0.89 30 0.82 0.77 0.84 0.79 Table 4: Summary of the two metrics NMI and OI for the different instance sizes. 5. Applications to Survey Data In the following section, we apply hypergraph modularity to find the hidden clusters of the public opinions emerging from survey data. In some social surveys, there are questions in which the respondents can elicit two or more items from a list. As discussed in Benati and Puerto (2024), the Eurobarometer is one example of these surveys. The Standard Eurobarometer includes questions about what are the two most important problems of the country. To that question, respondents can elicit at most two items from a closed list containing 17 terms, such as Terrorism, Unemployment, and so on. As described in Benati and Puerto (2024), the answers to that question can be modelled as a graph, in which nodes correspond to the items of the list, and there is an arc from item iand jfor every respondent that answered eliciting both iand j(note that multiple arcs and loops are allowed). The resulting graph, called the Items Graph, can be analysed using the network statistics and methodologies, such as modularity and community detection, and it has been found that network methods are much more informative than the alternative standard statistics, like k-means and principal component analysis. Indeed, modularity maximization found hidden opinion groups that could not be detected by k-means clustering or principal component analysis. In the following, we are applying modularity clustering to survey questions to which respondents can elicit more than two items. In the first example, coming from the Eurobarometer, respondents may elicit up to four items, chosen from a closed list. In the second example, respondents may answer up to three items, that can be chosen from a list and also freely pointed by respondents (forming an open list of items). We will show how hypergraph modularity can be applied to both settings and uncover hidden data patterns. A survey question in which respondents can elicit three (or more) number of items can be modelled as a hypergraph, in which nodes correspond to the list items, and there is a hyperarc between items i,jand kfor every answer eliciting the three items i, j, k. For example, we considered the Eurobarometer ZA7997, see European Commission (2024), held in June 2023. In that survey, there is a question in which respondents can elicit at most four items. The question is: On which of the following areas would you like the recovery plan NextGenerationEU to be spent in priority? The list proposed to respondents is: 14
•Scientific research and development •Youth, education and training •Culture and media •Energy •Transport •Climate change and environmental protection •Agriculture and rural development •Health •Improvement of working conditions of EU citizens •Support to SMEs •Defence and security •Public Administration •Industry •Digitalisation of economy and society •Other •None •Don’t know In a preliminary analysis, priorities could be grouped together into affinity areas. For example, we could expect that one group could consist of the economic priorities, such as Industry, Agriculture and Transport, while another group could be composed of the cultural investments, such as Research, Education and Culture. However, there are exponentially many possibilities by which the priorities can be a-priori combined and the cluster definitions can be biased by the researcher preferences. Therefore, some form of neutral, or non-arbitrary, data pre-processing is necessary: here, we apply the hypergraph modularity and model F2 z,γ. We consider two European nations, Spain and Italy, and we compare their expenditure priorities. When we analyse the frequencies by which priorities are mentioned, data are reported in Figure 6. It can be seen that there are some regularities as well as differences between the two nations. For example, Education is ranked high in both nations, while Health is an important priority to Spain, but not as much in Italy. Conversely, Energy, Defence and Digitalization are mentioned more often in Italy than in Spain. The Item Graph, see Benati and Puerto (2024), is a useful tool to visualize how answers are jointly selected. In Figure 7a, the hypergraph is projected into a simple graph: there is an arc between two items every time they were jointly selected by a respondent. In the original hypergraph, model (F2 z,γ) has been run with parameter γ=2 3: as respondents may elicit up to four items, a hyperarc is internal to a community under the following conditions. If the hyperarc contains two nodes, then both nodes must be in the same community. If it contains three nodes, then at least two of those nodes must be contained in the same community. If it contains four nodes, then at least three nodes must be contained in the same community. In Figure 7b, the item graph is drawn after the hypergraph clustering, clusters are composed of nodes with the same colour. As can be seen, the results of the two nations are different, but some regularity can be found. In both nations we can find that Health, Job condition, Environment, Education and Research are clustered together, forming the core of priorities that can be ascribed to redistributive policies, that we could term as Welfare spending priorities. It can be found, too, that Energy, Transport, SMEs, Defence, Public Administration and Industry are clustered together, forming the core of priorities that can be ascribed to Economic spending priorities. There are three items that are ambiguous: they are Agriculture, Culture, and Digitalization. The first two swapped their memberships between the two clusters, while Digitalization is a singleton cluster for Spain, but a Welfare issue for Italy. As shown in Benati and Puerto (2024), clusters can be used to reduce data dimension. Hyperedges correspond to respondents, so that we can analyse the feature of all the respondents within a cluster. In this example, respondents are classified as prioritizing the Welfare or the Economic spending, or termed as Unclassified if their hyperedges are not internal. After this reduction, in Figure 8a, it can be seen that in both countries the large majority of the citizens are in favour of the Welfare spending rather than the Economic spending, with the Spanish audience more oriented to it than the Italian. It is interesting to analyse further this finding. In political analysis, it is important to assess if political orientation plays a role on preferences. As can be seen in Figures 8b and 8c, it can be seen that in both countries the percentage of voters that are in favour of the Economic spending increases as the political preferences goes to the right of the political spectrum. In Spain, only 2% of left voters are in favour of the Economic spending, but they increase to 9% of the right voters. In Italy, they go from 8% of the left voters to 14% of the right voters. Conversely, the percentage of voters that in favour of the Welfare spending decreases. Therefore, hypermodularity clustering allowed to detect the regularity of the opinions along 15