A mixed-integer linear programming approach for soft graph clustering

TL;DR

This paper introduces a novel mixed-integer linear programming method for soft graph clustering that allows vertices to belong to multiple clusters with balanced memberships, overcoming limitations of previous approaches.

Contribution

It is the first approach to simultaneously allocate membership proportions and enforce balanced cluster memberships in soft graph clustering.

Findings

Clusters are not limited to k-clique neighborhoods.

Method produces non-trivial clusters in connected unweighted graphs.

Outperforms previous methods in flexibility and applicability.

Abstract

This paper proposes a Mixed-Integer Linear Programming approach for the Soft Graph Clustering Problem. This is the first method that simultaneously allocates membership proportion for vertices that lie in multiple clusters, and that enforces an equal balance of the cluster memberships. Compared to ([Palla et al., 2005], [Derenyi et al., 2005], [Adamcsek et al., 2006]), the clusters found in our method are not limited to k-clique neighbourhoods. Compared to ([Hope and Keller, 2013]), our method can produce non-trivial clusters even for a connected unweighted graph.

Tables

Table 1. Table 1: Computational results for KKI Instances with 20-50 vertices.

	Min Cut			Max Association
Data set	Opt	Gap	$r$	Opt	Gap	$r$
1019436-35	3.66	0.00	0.00	48.57	0.00	0.02
1541812-21	1.16	0.00	0.00	8.64	0.00	0.23
1577042-49	39.28	0.00	0.00	118.26	0.00	0.03
1638334-34	5.88	0.00	0.00	94.91	0.00	0.04
1735881-35	1.68	0.00	0.00	37.56	0.00	0.13
1779922-28	2.24	0.00	0.00	107.93	0.00	0.16
2371032-46	34.49	0.00	0.00	600.07	0.03	0.16
2558999-28	0.82	0.00	0.00	19.56	0.00	0.17
2601925-29	3.62	0.00	0.00	62.78	0.00	0.14
2618929-20	3.09	0.00	0.00	31.56	0.00	0.46
2768273-23	10.57	0.00	0.09	30.90	0.00	0.09
2903997-44	32.75	0.00	0.00	500.42	0.00	0.22
2917777-44	22.55	0.00	0.00	110.44	0.00	0.01
3310328-47	19.49	0.00	0.00	208.22	0.00	0.07
3611827-39	75.61	0.00	0.00	600.04	0.03	0.20
3713230-20	11.23	0.00	0.00	23.83	0.00	0.15
3902469-45	34.21	0.00	0.00	114.55	0.00	0.13
3972472-31	12.07	0.00	0.00	101.97	0.00	0.10
4104523-25	6.11	0.00	0.00	30.87	0.00	0.18
4275075-28	274.13	0.00	0.00	600.05	0.00	0.15
4362730-38	28.72	0.00	0.00	64.47	0.00	0.04
5216908-30	37.06	0.00	0.00	113.23	0.00	0.14
6346605-27	600.05	0.93	0.00	600.05	0.00	0.32
7129258-34	4.72	0.00	0.00	70.97	0.00	0.17
7415617-28	5.49	0.00	0.00	106.67	0.00	0.30
7774305-20	0.95	0.00	0.00	13.82	0.00	0.10
8263351-42	63.80	0.00	0.00	208.89	0.00	0.11

Table 2. Table 2: Computational results for problem instances with maximum edge weight M w = 50 subscript 𝑀 𝑤 50 M_{w}=50 and 100, K = 3 𝐾 3 K=3 clusters, and at least σ = 0.7 𝜎 0.7 \sigma=0.7 vertices must be in at least one cluster.

	Min Cut				Max Association
Data set	Opt	Gap	r	Con		Opt	Gap	r	Con
N15d015M50	0.939	0	0.215	100		2.347	0	0.769	100
(5/0)	0.269	0	0.195		(5/0)	0.987	0	0.155
N15d025M50	2.012	0	0.219	100		14.60	0	0.784	93.33
(5/0)	0.975	0	0.150		(5/0)	12.85	0	0.107
N15d05M50	78.46	0	0.393	100		412.3	0.074	0.912	100
(5/0)	43.20	0	0.317		(4/1)	85.48	0	0.198
N20d015M50	2.612	0	0.053	100		54.90	0	0.789	100
(5/0)	1.348	0	0.047		(5/0)	32.11	0	0.193
N20d025M50	66.07	0	0.146	100		252.2	0.067	0.876	100
(5/0)	51.50	-	0.110		(1/4)	0	0.017	0.177
N20d05M50	-	0.686	0.484	100		-	0.229	1.069	100
(0/5)	-	0.176	0.152		(0/5)	-	0.033	0.208
N30d015M50	31.38	0.826	0.025	100		-	0.044	0.817	100
(4/1)	11.07	0	0.048		(0/5)	-	0.021	0.073
N30d025M50	-	0.969	0.224	100		-	0.146	0.862	100
(0/5)	-	0.021	0.110		(0/5)	-	0.028	0.102
N50d015M50	-	1	0.142	93.33		-	0.176	0.937	100
(0/5)	-	-	0.095		(0/5)	-	0.032	0.066
N50d025M50	-	1	0.683	93.33		-	0.462	1.078	100
(0/5)	-	0	0.128		(0/5)	-	0.049	0.078
N15d015M100	0.540	0	0.089	80.00		3.193	0	0.484	80.00
(5/0)	0.492	0	0.123		(5/0)	5.069	0	0.476
N15d025M100	3.226	0	0.218	100		44.01	0	0.886	100
(5/0)	2.417	0	0.143		(5/0)	26.46	0	0.149
N15d05M100	35.86	0	0.359	100		342.5	0.116	1.070	100
(5/0)	20.11	0	0.202		(3/2)	95.53	0.021	0.184
N20d015M100	4.897	0	0.096	86.66		41.86	0	0.789	73.33
(5/0)	3.398	0	0.105		(5/0)	15.81	0	0.092
N20d025M100	71.12	0	0.102	100		261.9	0.049	0.922	100
(5/0)	62.69	0	0.076		(1/4)	0	0.018	0.146
N20d05M100	-	0.709	0.494	100		-	0.210	1.052	100
(0/5)	-	0.064	0.155		(0/5)	-	0.014	0.078
N30d015M100	46.06	0	0.006	93.33		515.5	0.037	0.731	93.33
(5/0)	70.07	0	0.007		(1/4)	0	0.024	0.123
N30d025M100	-	0.953	0.168	100		-	0.155	0.934	100
(0/5)	-	0.042	0.063		(0/5)	-	0.022	0.092
N50d015M100	-	1	0.129	93.33		-	0.187	0.905	100
(0/5)	-	0	0.089		(0/5)	-	0.099	0.095
N50d025M100	-	1	0.592	93.33		-	0.448	1.312	100
(0/5)	-	0	0.147		(0/5)	-	0.080	0.185

Table 3. Table 3: Computational results for problem instances with maximum edge weight M w = 50 subscript 𝑀 𝑤 50 M_{w}=50 , K = 2 𝐾 2 K=2 clusters, and at least σ = 0.7 𝜎 0.7 \sigma=0.7 vertices must be in at least one cluster.

	Min Cut				Max Association
Data set	Opt	Gap	r	Con		Opt	Gap	r	Con
N15d015M50	0.196	0	0.061	100		0.203	0	0.246	100
(5/0)	0.067	0	0.048		(5/0)	0.058	0	0.041
N15d025M50	0.180	0	0.014	90.00		0.256	0	0.264	100
(5/0)	0.038	0	0.027		(5/0)	0.066	0	0.079
N15d05M50	2.218	0	0.090	100		1.660	0	0.313	100
(5/0)	1.594	0	0.032		(5/0)	0.670	0	0.056
N20d015M50	0.211	0	0.014	90.00		0.699	0	0.214	100
(5/0)	0.067	0	0.032		(5/0)	0.288	0	0.051
N20d025M50	1.181	0	0.011	100		2.537	0	0.283	100
(5/0)	0.680	0	0.015		(5/0)	1.288	0	0.064
N20d05M50	46.63	0	0.122	100		61.40	0	0.352	100
(5/0)	29.13	0	0.064		(5/0)	26.19	0	0.037
N30d015M50	1.548	0	0.000	80.00		12.10	0	0.248	100
(5/0)	0.461	0	0.000		(5/0)	4.455	0	0.051
N30d025M50	63.92	0	0.023	100		169.4	0.006	0.323	100
(5/0)	71.98	0	0.015		(4/1)	185.0	0	0.050
N50d015M50	0	0.887	0.023	100		0	0.082	0.287	100
(0/5)	0	0.162	0.015		(0/5)	0	0.014	0.032
N50d025M50	0	0.997	0.111	100		0	0.177	0.354	100
(0/5)	0	0.002	0.028		(0/5)	0	0.014	0.055

Table 4. Table 4: Computational results for problem instances with maximum edge weight M w = 50 subscript 𝑀 𝑤 50 M_{w}=50 , K = 4 𝐾 4 K=4 clusters, and at least σ = 0.7 𝜎 0.7 \sigma=0.7 vertices must be in at least one cluster.

	Min Cut				Max Association
Data set	Opt	Gap	r	Con		Opt	Gap	r	Con
N15d015M50	3.453	0	0.566	75.00		92.08	0	1.614	65.00
(5/0)	2.826	0	0.343		(5/0)	68.72	0	0.260
N15d025M50	75.79	0	0.252	75.00		300.5	0.053	1.553	75.00
(5/0)	100.2	0	0.141		(2/3)	271.8	0.049	0.358
N15d05M50	0	0.617	0.863	75.00		0	0.232	1.737	75.00
(0/5)	0	0.090	0.294		(0/5)	0	0.015	0.236
N20d015M50	83.63	0	0.104	75.00		41.54	0.035	1.287	65.00
(5/0)	128.8	0	0.105		(1/4)	0	0.004	0.319
N20d025M50	0	0.724	0.322	75.00		0	0.136	1.821	75.00
(0/5)	0	0.252	0.169		(0/5)	0	0.020	0.143
N20d05M50	0	0.988	1.112	75.00		0	0.346	2.201	75.00
(0/5)	0	0.007	0.446		(0/5)	0	0.057	0.306
N30d015M50	237.4	0.943	0.060	75.00		0	0.080	1.670	75.00
(1/4)	0	0.085	0.070		(0/5)	0	0.029	0.180
N30d025M50	0	1	0.365	75.00		0	0.248	1.988	75.00
(0/5)	0	0	0.252		(0/5)	0	0.034	0.267
N50d015M50	0	1	0.794	65.00		0	0.447	1.955	70.00
(0/5)	0	0	0.468		(0/5)	0	0.102	0.162

Table 5. Table 5: Computational results for problem instances with maximum edge weight M w = 50 subscript 𝑀 𝑤 50 M_{w}=50 , K = 3 𝐾 3 K=3 clusters, and tested σ = 0.5 𝜎 0.5 \sigma=0.5 versus σ = 0.8 𝜎 0.8 \sigma=0.8 .

	Min Cut				Max Association
Data set	Opt	Gap	r	Con		Opt	Gap	r	Con
N15d015M50	0.574	0	0.075	100		0.698	0	0.363	100
(5/0)	0.305	0	0.100		(5/0)	0.206	0	0.334
N15d025M50	0.482	0	0.000	100		1.403	0	0.140	100
(5/0)	0.389	0	0.000		(5/0)	0.159	0	0.096
N15d05M50	24.29	0	0.258	100		116.2	0	0.391	100
(5/0)	12.43	0	0.174		(5/0)	33.96	0	0.116
N20d015M50	1.019	0	0.000	100		5.951	0	0.034	100
(5/0)	0.737	0	0.000		(5/0)	8.141	0	0.025
N20d025M50	12.31	0	0.018	100		149.4	0	0.237	100
(5/0)	10.26	0	0.016		(5/0)	199.2	0	0.083
N20d05M50	55.45	0.439	0.277	100		0	0.713	0.620	100
(1/4)	0	0.168	0.228		(0/5)	0	0.080	0.281
N30d015M50	3.175	0	0.000	93.33		357.0	0.813	0.086	93.33
(5/0)	2.334	0	0.000		(4/1)	196.5	0	0.067
N30d025M50	47.12	0	0.001	100		0	0.959	0.265	100
(5/0)	53.50	0	0.002		(0/5)	0	0.026	0.133
N50d015M50	209.7	1	0.007	93.33		0	1	0.313	100
(4/1)	212.0	0	0.016		(0/5)	0	0	0.137
N50d025M50	0	1	0.219	100		0	1	0.903	93.33
(0/5)	0	0	0.109		(0/5)	0	0	0.317

Equations

$$x_{i,c}$$

$$y_{i,c}$$

$$L_i$$

$$\sum _{c\in \mathcal{C}}{x}_{i,c}$$

$$(1-\delta )\sum _{i\in V}{x}_{i,c_1}\le \sum _{i\in V}{x}_{i,c_2}\le (1+\delta )\sum _{i\in V}{x}_{i,c_1},$$

$$\forall c_1,c_2\in \mathcal{C},c_1\ne c_2$$

$$y_{i,c_1}+y_{i,c_2}$$

$$t_{c_1,c_2}^i$$

$$\sum _{i\in V}{t}_{c_1,c_2}^i$$

$$t_{c_1,c_1}^i+t_{c_1,c_1}^j$$

$${\eta }_{c_1,c_2}^{i,j}$$

$$y_{i,c_1}+y_{j,c_2}+a_{ij}+(1-{\eta }_{c_1,c_2}^{i,j})$$

$$s_{c_1,c_2}^e$$

$$v_{c_1,c_2}^e=\left(x_{i,c_1}+x_{j,c_2}\right)s_{c_1,c_2}^e$$

$${\tau }_{c_1,c_2}^{e,i}$$

$${\tau }_{c_1,c_2}^{e,j}$$

$$\kappa (c_1,c_2)=\sum _{(i,j)\in E}{w}_{i,j}({\tau }_{c_1,c_2}^{e,i}+{\tau }_{c_1,c_2}^{e,j})$$

$$y_{i,c}+y_{j,c}+a_{ij}$$

$$z_{i,j}^c$$

$$A(c)=\sum _{(i,j)\in E}{w}_{i,j}(x_{i,c}+x_{j,c})z_{i,j}^c$$

$${\pi }_{i,j}^{c,i}$$

$${\pi }_{i,j}^{c,j}$$

$$A(c)=\sum _{(i,j)\in E}{w}_{i,j}({\pi }_{i,j}^{c,i}+{\pi }_{i,j}^{c,j})$$

$${\gamma }_{i,j}^c$$

$$\sum _{i\in V}{y}_{i,c}-1$$

Full text

1 Introduction

Consider an undirected graph $G=(V,E)$ with $V$ the set of vertices and $E=\{(i,j)\ |\ i,j\in V,i<j\}$ the set of edges. Each edge $e\in E$ is associated with a weight $w_{e}\in\mathbb{R}$ that indicates the similarity between its two end vertices–the lager the weight, the more “similar” the two vertices are. The hard graph clustering (HGC) problem is to create distinct partitions (clusters, or, communities) of the set of vertices according to their similarities, i.e., to form $V_{1},\ldots,V_{k}$, where $\bigcup_{i=1,\ldots,k}V_{i}=V$, and $V_{i}\cap V_{j}=\emptyset$ for all $i,j\in\{1,\ldots,k\}$, $i\neq j$. For a thorough literature review of the graph clustering problems, see, e.g., Schaeffer (2007), and for fast algorithms for large-scale networks, see, e.g., Girvan and Newman (2002); Clauset et al. (2004); Rosvall and Bergstrom (2008). For datasets, see, e.g., SNAP network datasets sna and SNAP biomedical datasets Marinka Zitnik and Leskovec (2018).

The soft graph clustering (SGC) problem, (also known as fuzzy graph clustering), on the other hand, allows clusters to have overlaps. A vertex may be a member of more than one cluster. There are numerous applications of SGC, such as: brain research, social network research, natural language processing, citation, and collaboration networks, and so on. A precise problem definition of the SGC varies and is dependent on the application, and sometimes it may not be possible to provide a precise problem definition.

The subject of study in this paper considers the combinatorial optimisation problem where we are required to determine: 1) the composition of each of the clusters; 2) for each vertex that belongs to more than one cluster, how the membership is distributed amongst the clusters (we denote this by $x_{ic}$, for Vertex $i$ in Cluster $c$, hence $\sum_{c}x_{ic}=1$ for all $i\in V$). We consider the case that an equal balance of the cluster total vertex memberships is desirable, and that not all vertices are required to be in a cluster. We consider two equally important objectives: 1) to minimize the sum of inter-cluster edge weights (cut across clusters); and 2) to maximize the sum of intra-cluster edge weights (cluster association).

1.1 Literature review

There are a number of existing soft clustering algorithms, each designed to suit different applications see, e.g., CFinder of Palla et al. (2005), (see also Derényi et al. (2005), Adamcsek et al. (2006)), the MaxMax Algorithm of Hope and Keller (2013), the WATSET methods of Ustalov et al. (2018) for NLP, the Chinese Whisper method of Biemann (2006), Betweenness-based method of Pinney and Westhead (2006), and the Purifying and Filtering the Coupling Matrix approach of Liu and Foroushani (2016). Of these methods, Biemann (2006), Pinney and Westhead (2006), Liu and Foroushani (2016), and Ustalov et al. (2018) are designed for unweighted graphs only (i.e., graphs with unit edge weight).The MaxMax Algorithm is designed for weighted undirected graphs. For unweighted graphs, however, it will return a trivial solution–each connected component of the graph will be a cluster. The CFinder is based on the finding of $k$-clique neighbourhoods. The mixed-integer linear programming method we propose in this paper is able to accomodate both weighted and unweighted graphs, with a small modification required for the latter. In the preliminary experiments section, we will compare and contrast the different methods.

We are not aware of any mixed-integer linear programming (MILP) models for the SGC problems. There are, however, MILP models for other graph clustering, machine learning, and data classification problems. The article Bertsimas and Shioda (2007) presents MILP formulations for classification and regression. The idea for classification, e.g., is to partition Class 1 points into $K$ disjoint subsets by finding the hyperplanes that describe the partitioning polyhedrons such that no Class 0 points can be expressed as a convex combination of the Class 1 points in each partition. In general clustering problems, Sağlam et al. (2006) proposes a MILP formulation where one wishes to partition a set of data set into $k$ (a predetermined number of) clusters. The objective is to minimize the maximum diameter of the generated clusters in order to obtain evenly compact clusters. Essentially the method is an IP-based heuristic method, with some variables fixed by the solution of maximal independent set of size $k$, where each member of this set is a seed member in the $k$ clusters. The IP model (which is in fact a bilinear model, but linearized using standard linearization strategies) is then solved to obtain an optimal solution to the general clustering problem. Other MILP-based work can be found in, e.g., Gilpin et al. (2013) and Ye (2007) for hierarchical clustering. The latter presents an application in recommendation systems. In Clique Covering Problem (CCP), an NP-hard combinatorial optimisation problem where an undirected graph is to be partitioned to form complete subgraphs, Miyauchi et al. (2018) proposes a compact ILP formulation for a relaxed problem, as well as a post-optimization repair procedure and a proof of optimality for the final solution to the original problem.

1.2 Contribution of the paper

As far as we aware, this paper is the first to propose a methodology for SGC that i) deals with undirected graphs with general integer edge weights $w_{e}\geq 0$ and truely takes the values of $w_{e}$ into the optimization process in the way that the larger the value is, the more favourably it will be considered, and at the same time, also deals with unweighted graphs (i.e., graphs with unit edge weight); ii) simultaneously allocates membership proportion (the $x_{ic}$ values) for vertices that lie in multiple clusters; and iii) it enforces an equal balance of the clusters that the sum of vertex memberships over all clusters are roughly the same. We propose an approach that is based on a polynomial-size MILP model, beginning with a small value for $K$–we enforce that the graph has at least $K$ clusters. One can apply an adaptive approach to find the best value of $K$ iteratively, but the focus of this research is to solve an instantaneous SGC problem with a given $K$.

The method of Palla et al. (2005) requires obtaining $E^{\prime}=\{e\in E\ |\ w_{e}>w^{*}\}$, and the graph is subsequently clustered by finding $\kappa$-clique neighbourhoods on an unweighted graph $H=(V,E^{\prime})$. Our method does not require the finding of $\kappa$-cliques, and it takes the values of $w_{e}$ into account during optimization. Comparing with the method of Hope and Keller (2013) (MaxMax), for unweighted graphs, if the graph is connected, then MaxMax will produce only one cluster which is the entire graph. Our method, however, can deal with unweighted graphs by converting them into weighted ones via a simple transformation.

In Section 2.1, we present our basic MILP model by considering a number of standard requirements for the SGCP. In Section 2.2, we discuss our strategies for graph connectivity. In Section 2.3, we discuss two objectives: 1) minimizing the total inter-cluster cut, and 2) maximize the total intra-cluster association. In Section 3, we present preliminary numerical results. We then conclude our findings and discuss future research directions in Section 4.

2 A mixed-integer linear programming formulation

2.1 The basic model

We first introduce the notation used in this paper. Let:

• •

$A=\{a_{ij}\in\{0,1\}\ |\ i,j\in V\}$ be the adjacency matrix of $G$;

• •

$M_{w}=\max\{w_{e}\ |\ e\in E\}$ the maximum edge weight;

• •

$K$ be the number of clusters;

• •

$\mathcal{C}=\{c_{1}\ldots,c_{K}\}$ be the set of clusters;

• •

$y_{i,c}\in\{0,1\}$ be a binary decision variable with $y_{i,c}=1$ indicating Vertex $i$ is a member of Cluster $c$;

• •

$x_{i,c}\in[0,1]$ be a continuous decision variable indicating the membership of Vertex $i$ in Cluster $c$;

• •

$\kappa(c_{1},c_{2})$ be the cut between clusters $c_{1}$ and $c_{2}$;

• •

$0<\mu<1$ a predetermined minimum membership if a vertex is a member of a cluster;

• •

$0<\delta<1$ a predetermined tolerance equal balance of cluster membership;

• •

$0<\nu<1$ a predetermined maximum overlap factor.

Now we introduce the constraints. First, we have a set of Membership Constraints. The membership of Vertex $i$ in Cluster $c$ can only be non-zero if it is a member of $c$, and when it is, the membership must be no less than a predetermined value.

[TABLE]

Let $L_{i}$, for each $i\in V$, be an auxiliary binary variable with $L_{i}=1$ indicating $i$ is in at least one of the clusters and $L_{i}=0$ otherwise. We require that the sum of memberships for any vertex over all clusters is exactly 1 if the vertex is a member of at least one cluster and 0 otherwise.

[TABLE]

We consider an Equal Balance Requirement where the sum of memberships in all clusters are “roughly” the same. As far as we aware, this is the first method that considers such a requirement. This requirement can be modelled as below.

[TABLE]

Next, we consider the Overlap Cardinality Constraints. Let $t^{i}_{c_{1},c_{2}}$ be auxiliary binary decision variable such that $t^{i}_{c_{1},c_{2}}=1$ if and only if $i$ is a member of both clusters $c_{1}$, $c_{2}$, i.e., $y_{i,c_{1}}=y_{i,c_{2}}=1$. We have that:

[TABLE]

Given any pairs of clusters $c_{1},c_{2}\in\mathcal{C}$, with $c_{1}\neq c_{2}$, the number of vertices that are in both clusters cannot be larger than a predetermined fraction, $\nu$, of the cardinality of either of the clusters.

[TABLE]

To calculate the Inter-cluster Cuts between $c_{1}$ and $c_{2}$, we require auxiliary binary variables and nonlinear terms. First, let $\eta^{i,j}_{c_{1},c_{2}}$ be a binary variable such that $\eta^{i,j}_{c_{1},c_{2}}=1$ if both of $i,j$ are in the intersection of clusters $c_{1}$ and $c_{2}$.

[TABLE]

We then use a binary variable $s^{e}_{c_{1},c_{2}}$ to indicate the existence of an edge (cut) $e=(i,j)$ (i.e., $a_{ij}=1$) with $i$ in $c_{1}$ and $j$ in $c_{2}$, but not both in the intersection of $c_{1}$ and $c_{2}$ (otherwise the “cut” should not be counted). The constraints are as below. For each pair of distinct vertices $i,j\in V,\ i\neq j$ and each pair of distinct clusters $c_{1},c_{2}\in\mathcal{C},\ c_{1}\neq c_{2}$, we have that:

[TABLE]

Now, the cut between two distinct vertices $i,j\in V$ across two distinct clusters $c_{1},c_{2}\in\mathcal{C}$, if edge $e=(i,j)\in E$ exists, is defined by: $w_{e}v^{e}_{c_{1},c_{2}}$, for

[TABLE]

The terms $x_{i,c_{1}}s^{e}_{c_{1},c_{2}}$ and $x_{j,c_{2}}s^{e}_{c_{1},c_{2}}$ are bilinear, and can be linearized by introducing auxiliary non-negative continuous variables $\tau^{e,i}_{c_{1},c_{2}}\geq 0$ and $\tau^{e,j}_{c_{1},c_{2}}\geq 0$ and the following constraints. For each $e=(i,j)\in E,\ c_{1},c_{2}\in\mathcal{C},\ c_{1}\neq c_{2}$, we have that:

[TABLE]

The cut $\kappa(c_{1},c_{2})$ is calculated by the following linear term.

[TABLE]

Now, we consider the Intra-cluster Association calculations. Let $z^{c}_{i,j}$ be an auxiliary binary variable with $z^{c}_{i,j}=1$ if $i,j$ are both in $c$ and that the edge $e=(i,j)$ exists in $E$. The constraints are give by:

[TABLE]

The intra-cluster association of $c\in\mathcal{C}$ is given by

[TABLE]

We linearize the association using auxiliary continuous variables $\pi^{c,i}_{i,j}$ and $\pi^{c,j}_{i,j}$ to capture the memberships of two vertices in the same cluster should an edge exists between them. I.e., when $z^{c}_{i,j}=1$, $\pi^{c,i}_{i,j}=x_{i,c}$, otherwise, $\pi^{c,i}_{i,j}=0$. (Similarly for $\pi^{c,j}_{i,j}$). The constraints are as below.

[TABLE]

The association within a cluster $c$, denoted by $A(c)$, is calculated as follows:

[TABLE]

2.2 Cluster connectivity

Clusters are required to be connected, it is likely that this can only be achieved with exponentially many variables and solved using a branch-and-price approach, or with exponentially many constraints and solved using a branch-and-cut approach. Our approach here is not an exact one, as we wish to keep the formulation compact (i.e., polynomial in size). We derive a number of constraints that capture a few conditions for graph connectivity that are necessary but not sufficient. When there is a connectivity violation, violation elimination constraints can be added in a lazy fashion. (Lazy constraints is a technical term in integer programming–hard (and often exponentially many) constraints are relaxed, and are only added when the current integer optimal solution to the relaxed problem violates them–usually only a very small number of them are violated–the problem is then re-optimized, and the procedure recurs until there is no more violated hard constraints). In our preliminary test where we used randomly generated undirected graphs (see Preliminary Numerical Results section), most of the problem instances produced connected clusters.

First of all, the number of edges in a connected undirected graph cannot be smaller than the cardinality of the graph minus one. We define $\gamma^{c}_{i,j}\in\{0,1\}$ to be a span variable that can only be one when $i$ and $j$ are both in Cluster $c$ and that the edge $(i,j)$ exists in $E$. We have that:

[TABLE]

and that

[TABLE]

We also require that all vertices in each cluster must be connected to at least one other vertex in the same cluster.

[TABLE]

However, Constraints (39)–(43) are not enough to eliminate multiple loops within a cluster. A cluster may contain vertices $\{1,2,3,4,5,6\}$, but the constraints cannot eliminate the formation of subgraphs $\{1,2,3\}$ and $\{4,5,6\}$ within the cluster.

Therefore, we borrowed the time constraints idea for Asymmetric Travelling Salesman Problem (ATSP) subtour elimination (Miller et al. (1960)). Let $t_{i}\geq 0$, $i\in V$, be a decision variable indicating the “time of arrival” at Vertex $i$. Suppose $\{1,2,3\}$ is a strict subset of a cluster $c$, to prevent the relevant span variables from forming a loop, we require that $\gamma^{c}_{i,j}=1$ if and only if $t_{j}=t_{i}+1$, for each pair of distinct $i,j\in V$, $i\neq j$ and each $c\in\mathcal{C}$. We have that:

[TABLE]

E.g., consider the vertices $\{1,2,3\}$ a strict subset of $c$, a loop with edges $(1,2),(2,3),(1,3)$ a violation of the time constraints, as $t_{3}$ cannot be equal to $t_{1}+1$ and $t_{1}+2$ simultaneously. Constraints (44) and (45) can eliminate certain types of loops formed by the span variables, but not all. In any case, adding them will give the span variables a better chance at making a full span in a cluster and thus the connectivity.

Notice that the time constraints do not cut off feasible solutions due to the fact that even if a variable $\gamma^{c}_{i,j}$ is forced to be 0, there is no impact on the values of $y_{ic}$ and $y_{jc}$, as (39)–(41) do not induce a bi-conditional relation. In our preliminary numerical experiments, the time constraints are expensive to implement, and often we do not have disconnected clusters with just (39)–(43) alone.

When cluster size $K$ is a hard constraint, and there exist a disconnected cluster, in the case of maximizing total association, e.g., one can consider adding the following cut in a lazy fashion and re-optimize. Let $\mathcal{I}^{*}$ be the set of $y$-variables with a value of 1 in the optimal solution, we add:

[TABLE]

2.3 The objective function

Two commonly considered objective functions are: 1) minimize the sum of inter-cluster cuts and 2) maximize the sum of intra-cluster association given by (47) and (48) below respectively.

[TABLE]

[TABLE]

Notice that it is necessary to enforce a minimum cluster size, otherwise a trivial optimal solution to (47) is to assign no vertices to any clusters. Thus we have:

[TABLE]

for $0<\sigma<1$ a predetermined value.

The two objectives should be considered simultaneously. Some may consider minimizing the sum of ratio of inter-cluster cuts and intra-cluster associations over all pairs of distinct clusters.

[TABLE]

Diagram (1) below shows how the sum of ratios of inter-cluster cuts and intra-cluster associations changes. The first data point on the left is obtained when we optimize (47) and set $w^{1}$ to be the value of $\sum_{\forall c\in\mathcal{C}}A(c)$ in the optimal solution. We can see that the total intra-cluster association is also small, and so is the value of (50). The last data point on the right is obtained when we optimize (48) and we set $w^{2}$ to be the optimal objective value. We then minimize (47) with $\sum_{\forall c\in\mathcal{C}}A(c)\geq\ell_{j}$, for $\ell_{j}=\frac{j}{10}(w^{2}-w^{1})$, $j=1,\ldots,10$.

One can, however, consider recent advances in Bi-objective Integer Programming (see, e.g., Dai and Charkhgard (2018)). In fact, different applications of the SGC may have a different objective function that describes the SGC better. Besides, the structure of the graph must be taken into consideration in determining what the most appropriate objective function is. One may therefore consider applying machine learning to automate the finding of an appropriate objective function.

2.4 Edge weight transformation for unweighted graphs

For undirected graphs with $w_{e}=1$, for all $e=(i,j)\in E$, the MILP does not work very well because it cannot distinguish between an edge in a sparse neighbourhood with one in a dense neighbourhood. To give favour to edges in a dense neighbourhood, we obtain new edge weights by calculating $w^{\prime}_{e}=1+|\{k\in V\ :\ (i,k),(j,k)\in E,\ k\neq i,j\}|$, i.e., the new edge weight of $e$ will be one plus the number of vertices that are connected to both of the end nodes of $e$.

3 Preliminary numerical results

We compared our method with CFinder Palla et al. (2005); Derényi et al. (2005); Adamcsek et al. (2006) and MaxMax Algorithm Hope and Keller (2013). We used a 21-vertex instance of weighted graph. We can see that the clusters of CFinder are the $3$-clique neighbourhoods, and therefore some of the edges with heavy weight are not included, e.g., $(6,7)$, $(9,10)$, $(17,20)$.

The MaxMax solution, however, have these heavy-weight edges covered. However, for unweighted graphs, MaxMax will cluster the graph by connected components, so if the graph is connected, then there will be only one cluster.

The MILP solution produces not only the clusters, but also the membership proportion for each vertex that belongs to more than one cluster. None of MaxMax or CFinder provides this information. One can see clearly that the Min Cut solution can also produce a relatively balanced set of clusters in terms of total weighted memberships of the clusters when the clusters are connected. None of the existing SGC methods considered such a constraint.

In the Max Association solution below, unfortunately, because the MILP enforced three clusters, and the blue cluster is not connected, we have four clusters instead, so the equal balance cannot be guarantee in this case. In any case, the four clusters demonstrated the four strongest connections by $w_{e}$. As for the membership proportions, take Vertex 9 as an example, since there are more heavily weighted links in the pink cluster rather than the green cluster, 0.9 of its membership is in the former, and 0.1 in the latter.

3.1 KKI instances

In this section, we applied our method on some KKI instances from https://github.com/shiruipan/graph_datasets. Since the KKI instances are unweighted, we perform the edge weight transformation described in Section 2.4.

In Table 1, column Opt presents the time taken (in seconds) for proving optimality, and the column Gap presents the gap between the MIP relative gap, (i.e., the gap between the best integer objective and the objective of the best active node on the branch-and-bound tree), for instances that are not solved to optimality within the given time limit (10 minutes). The column r is the ratio of total inter-cluster cut over total intra-cluster association. We used IBM ILOG CPLEX Optimization Studio v12.8 Support (2018) for solving the MILPs.

3.2 Instances from a random problem generator

We have generated 10 problem classes with 5 instances for each problem class. The problem instances were generated using a random problem generator where number of vertices on the undirected graph, ($|V|$), density of graph $\frac{|E|}{2|V|(|V|-1)}$, and an upper bound on the edge weight $M_{w}$ are taken as user input. The edge weights are generated by each given a number chosen uniformly randomly between $[1,M_{w}]$. In Table 2, the names of the problem classes are in the form of: N followed by the value of $|V|$, d followed by the density of $G$, and then $M$ followed by the value of $M_{w}$. The column CON is the percentage of clusters that are connected in each problem class. In the columns where two numbers were given for each problem class, the first row presents the average value over 5 problem instances, and the second row presents the value of the standard deviation. The bracket underneath each problem class indicates the number of instances that are solved to optimality versus the number of instances that are not.

3.3 Preliminary testing

We experimented with instances with maximum edge weight $M_{w}=50$ and $100$, $K=3$ clusters, and $\sigma=0.7$. From the results of Table 2, we can see that the computation time grows exponentially as the size of the problem grows, or the density of the graph grows. For problem instances that are not solved to optimality, the MIP gaps are large. However in all problem instances, feasible solutions are found reasonably quickly. We can see that the computation time for maximizing total intra-cluster association is substantially longer than minimizing total inter-cluster cuts. For problem instances that are not solved to optimality, however, the former has a smaller MIP gap whilst the latter can solve larger problem instances. With an objective of maximizing total association, we did not include Constraint (49) as the objective will drive as much edges used as possible. We cannot obtain a conclusive remark on the effects of $M_{w}$ in terms of computation time.

4 Conclusions and future research directions

We proposed and developed a polynomial-size MILP model for the SGCP. As far as we aware, this is the first approach that simultaneously allocates membership proportion (the $x_{ic}$ values) for vertices that lie in multiple clusters, and that enforces an equal balance of the clusters so that the sum of vertex memberships over all clusters are roughly the same. Compared to CFinder (Palla et al. (2005), Derényi et al. (2005), Adamcsek et al. (2006)), the clusters found in our method are not limited to $k$-clique neighbourhoods. Compared to MaxMax (Hope and Keller (2013)), our method can produce non-trivial clusters even for a connected unweighted graph. An obvious future research direction is to perform a thorough numerical experiment on all parameters used in the model. One can consider alternative formulations, e.g., a branch-and-cut approach, using constraints to cut off infeasible solutions (namely, solutions with disconnected clusters); or, a branch-and-price approach, using exponentially many variables each representing a feasible cluster (in which case connectivity is guaranteed). Even though the cut separation and column generation subproblems themselves are expected to be NP-hard, heuristic approaches can be derived for speedy execution. Further, one may consider an adaptive approach to find the optimal number of clusters, instead of iteratively optimising over different values of $K$.