Mixed-integer programming techniques for the minimum sum-of-squares clustering problem

Burgard, Jan Pablo; Moreira Costa, Carina; Hojny, Christopher; Kleinert, Thomas; Schmidt, Martin Article — Published Version Mixed-integer programming techniques for the minimum sum-of-squares clustering problem Journal of Global Optimization Provided in Cooperation with: Springer Nature Suggested Citation: Burgard, Jan Pablo; Moreira Costa, Carina; Hojny, Christopher; Kleinert, Thomas; Schmidt, Martin (2023) : Mixed-integer programming techniques for the minimum sum-of-squares clustering problem, Journal of Global Optimization, ISSN 1573-2916, Springer US, New York, NY, Vol. 87, Iss. 1, pp. 133-189, https://doi.org/10.1007/s10898-022-01267-4 This Version is available at: https://hdl.handle.net/10419/307510 Standard-Nutzungsbedingungen: Die Dokumente auf EconStor dürfen zu eigenen wissenschaftlichen Zwecken und zum Privatgebrauch gespeichert und kopiert werden. Sie dürfen die Dokumente nicht für öffentliche oder kommerzielle Zwecke vervielfältigen, öffentlich ausstellen, öffentlich zugänglich machen, vertreiben oder anderweitig nutzen. Sofern die Verfasser die Dokumente unter Open-Content-Lizenzen (insbesondere CC-Lizenzen) zur Verfügung gestellt haben sollten, gelten abweichend von diesen Nutzungsbedingungen die in der dort genannten Lizenz gewährten Nutzungsrechte. Terms of use: Documents in EconStor may be saved and copied for your personal and scholarly purposes. You are not to copy documents for public or commercial purposes, to exhibit the documents publicly, to make them publicly available on the internet, or to distribute or otherwise use the documents in public. If the documents have been made available under an Open Content Licence (especially Creative Commons Licences), you may exercise further usage rights as specified in the indicated licence. https://creativecommons.org/licenses/by/4.0/ Journal of Global Optimization (2023) 87:133–189 https://doi.org/10.1007/s10898-022-01267-4 Mixed-integer programming techniques for the minimum sum-of-squares clustering problem Jan Pablo Burgard1·Carina Moreira Costa2·Christopher Hojny3· Thomas Kleinert4·Martin Schmidt2 Received: 9 March 2022 / Accepted: 13 December 2022 / Published online: 10 January 2023 © The Author(s) 2023, corrected publication 2023 Abstract The minimum sum-of-squares clustering problem is a very important problem in data mining and machine learning with very many applications in, e.g., medicine or social sciences. However, it is known to be NP-hard in all relevant cases and to be notoriously hard to be solved to global optimality in practice. In this paper, we develop and test different tailored mixed-integer programming techniques to improve the performance of state-of-the-art MINLP solvers when applied to the problem—among them are cutting planes, propagation techniques, branching rules, or primal heuristics. Our extensive numerical study shows that our techniques significantly improve the performance of the open-source MINLP solver SCIP. Consequently, using our novel techniques, we can solve many instances that are not solvable with SCIP without our techniques and we obtain much smaller gaps for those instances that can still not be solved to global optimality. Keywords Minimum sum-of-squares clustering ·Mixed-integer nonlinear optimization · Global optimization ·Computational techniques BMartin Schmidt [email protected] Jan Pablo Burgard [email protected] Carina Moreira Costa [email protected] Christopher Hojny [email protected] Thomas Kleinert [email protected] 1Department of Economic and Social Statistics, Trier University, Universitätsring 15, 54296 Trier, Germany 2Department of Mathematics, Trier University, Universitätsring 15, 54296 Trier, Germany 3Combinatorial Optimization Group, Department of Mathematics and Computer Science, Eindhoven University of Technology, P.O. Box 513, 5600 MB Eindhoven, The Netherlands 4Quantagonia GmbH, Foellerweg 37, 61352 Bad Homburg, Germany 123 134 Journal of Global Optimization (2023) 87:133–189 Mathematics Subject Classification 90C10 ·90C11 ·90C57 ·90-08 1 Introduction Given a set of data points in a normed vector space and a number of clusters, the clustering problem consists in deciding which data point should be assigned to which cluster. Moreover, a representative point for each cluster needs to be determined. Clustering problems form a highly relevant sub-class of unsupervised learning in machine learning and computational statistics. Its relevance is supported by many applications, e.g., in functional data analysis [10,53], image processing [9], bio-informatics [12], economics [32], and social sciences [31]. For a detailed survey of the history of clustering problems we refer to Steinley [58]. Depending on, e.g., the way how vicinity is measured and on whether the representative is an arbitrary point or one of the data points, different variants of clustering problems arise. In this paper we consider the minimum sum-of-squares clustering (MSSC) problem. Here, distance between data points is measured using the squared Euclidean norm and any point can be chosen as the representative for each cluster. Modeling this problem leads to a nonconvex mixed-integer nonlinear optimization problem (MINLP) that is extremely hard to solve for high-dimensional real-world problems. Moreover, the problem is known to be NP-hard even in the case of two dimensions; see, e.g., Aloise et al. [11], Dasgupta [2], and Mahajan et al. [41]. This is why the problem is most frequently solved using heuristics out of which the k-means clustering method is the most prominent one; see, e.g., Lloyd [39], MacQueen [40]. However, solving such clustering problems only heuristically may come with severe disadvantages. Since solving the MSSC problem is an unsupervised learning problem, the outcome typically requires the interpretation of experts from the specific field of application such as medicine or social sciences. This interpretation, however, may be completely wrong in the case that the expert is confronted with a heuristic clustering solution of bad quality. Moreover, it is easy to imagine that such a misleading interpretation might have some severe, e.g., medical, consequences. Thus, there is a strong need for sophisticated optimization techniques to improve the process of solving clustering problems to global optimality and this is exactly the contribution of this paper: We take the MINLP solver SCIP and enhance its solution process by developing novel mixedinteger optimization techniques that enable us to solve MSSC problems to global optimality that cannot be solved with the plain version of SCIP. Of course, we are not the first ones trying to solve the MSSC problem to global optimality. To the best of our knowledge, the earliest application of branch-and-bound methods is presented by Fukunaga et al. [25], which has been refined later on by Diehr [15]. A variant of a so-called repetitive branch-and-bound method has been devised by Brusco [7], where the authors conclude that their method is well-suited for a small number of clusters. Another recent branch-and-bound approach is presented by Sherali and Desai [55]. The authors use reformulation-linearization-techniques (RLT) embedded in a branch-and-bound method to solve the problem to global optimality. In their introduction, they also talk about a “limited number of optimization techniques” as opposed to a rather large number of heuristics that are used in practice. As an additional technique, the authors present further valid inequalities to tackle the inherent symmetry of the problem. Regarding symmetry breaking for clustering problems we also refer to Plastria [48], which is, generally speaking, a modeling tutorial paper but which also contains a discussion of symmetry breaking constraints for the cluster123 Journal of Global Optimization (2023) 87:133–189 135 ing problem. Aloise and Hansen [4] also consider the MSSC problem and try to re-produce the results of Sherali and Desai [55]. However, the re-production failed since significantly longer running times have been observed. Consequently, the computational efficiency of the RLT-based branch-and-bound method may be taken with some care. Another algorithmic technique is the column generation approach presented first in Merle et al. [18] and which has been re-considered and improved by Aloise et al. [5]. Also other classic techniques of mixed-integer (non)linear optimization have been applied such as generalized Benders decomposition in Floudas et al. [21], and Tan et al. [59]. Alternatively, Peng and Xia [45] consider the MSSC problem as a concave minimization problem and adapt Tuy’s cut method, see Horst and Tuy [34], for solving the problem. Further, Prasad and Hanasusanto [49]propose improved conic reformulations of the MSSC problem and also study some symmetry breaking techniques. Tîrn˘auc˘aetal.[60] follow a more geometric approach that is based on Voronoi diagrams. Finally, there is a rather large branch of literature towards the application of techniques from semi-definite programming (SDP). Peng and Wei [46], Peng and Xia [44] proved the equivalence between the MSSC problem and a 0-1 SDP reformulation. Based on this 0-1 SDP model, Aloise and Hansen [3] propose a branch-and-cut algorithm and solve instances with up to 202 data points to global optimality. More recently, Piccialli et al. [47] consider the same mixed-integer SDP for the MSSC problem and propose another branch-and-bound algorithm that is capable of solving real-world instances with up to 4000 data points. To the best of our knowledge, this is the most recent state-of-the-art branch-andbound algorithm for the MSSC problem. SDP-like models have also been used in De Rosa and Khajavirad [13]touseZ=XX∈[0,1]n×nfor encoding the clustering instead of X∈{0,1}n×k. The authors derive cutting planes and show a relation of their cutting planes to the cut polytope; see Deza and Laurent [14] for a survey on the latter. The presented numerical experiments show that these novel cutting planes can be strong but the authors only solve the initial LP relaxation and do not apply a complete branch-and-bound method. Finally, some recent ideas based on reduced-space techniques seem to be very promising, see Hua et al. [35] and Liberti and Manca [38]. Besides that, in Liberti and Manca [38], the authors discuss the MSSC problem with several side constraints. One of their base models is, in particular, the convex MINLP that we present in the next section. In our contribution, we add to the literature on solving the MSSC problem to global optimality. To this end, we develop novel mixed-integer programming techniques that are mainly motivated by geometric insights and that improve the branch-and-cut solution process of an MINLP solver. To be more precise, we present two MINLP formulations of the problem (Sect.2), develop cutting planes (Sect.3), propagation methods (Sect.4), as well as problemspecific branching rules (Sect.5) and primal heuristics (Sect.6). We implement and test all techniques in the open-source MINLP solver SCIP; see Gamrath et al. [26]. By doing so, we also automatically apply state-of-the-art symmetry breaking techniques to the problem. Our numerical results are presented and discussed in Sect.8, where we show that our techniques significantly improve the solution process. We close the paper with some concluding remarks and some potential topics for future work in Sect.9. Our code is publicly available at GitHub.1 Although our numerical results clearly show that the solution process of an MINLP solver applied to the MSSC problem is significantly improved, we do not beat current state-ofthe-art and SDP-based techniques as studied in Piccialli et al. [47]. Nevertheless, we are convinced that it is worth to also push MINLP-based approaches forward so that, in the end, different techniques for various approaches can be combined to lead to an even better and maybe hybrid solution approach. 1https://github.com/christopherhojny/globally-solving-MSSC. 123 136 Journal of Global Optimization (2023) 87:133–189 2 MINLP models for the MSSC problem We now model the minimum-sum-of-squares clustering (MSSC) problem as a mixed-integer nonlinear optimization problem (MINLP). To this end, we are given a set of data points p∈P⊆Rdand a positive integer 2 ≤k≤|P|, which is the number of clusters of the problem. The task is then to assign every data point p∈Pto a cluster (indexed by j∈[k]:={1,...,k}) so that the sum of the squared Euclidean distances between the data points and the corresponding centroids cjis minimal. This problem is modeled via the following MINLP: min x,c p∈P j∈[k] xpjp−cj2(1a) s.t. j∈[k] xpj =1,p∈P,(1b) xpj ∈{0,1},p∈P,j∈[k],(1c) cj∈B,j∈[k].(1d) The binary variables xpj are the assignment variables that model whether the data point pis assigned to cluster j(xpj =1)or not (xpj =0). Moreover, B⊆Rdis a set that contains all points in P. This can, e.g., be the bounding box of P. That is, if for each i∈[d], i=min{pi:p∈P}and ui=max{pi:p∈P},thenB={c∈Rd:i≤ci≤ ui,i∈[d]} is a valid choice. Note that (1d) is not necessary for the correctness of Model (1). Nevertheless, we include it in our implementation, because Model (1) is a nonconvex MINLP, for which bounds on variables are usually beneficial. The objective function measures the sum of the squared Euclidean distances between the data points and the centroids of the clusters to which they belong. Finally, Constraint (1b) ensures that every point is assigned to exactly one cluster. Note that this model is cubic since the objective function uses multiplications of the assignment variables xwith the norms that depend on the centroids c, which are variables of the problem as well. In particular, Model (1) is a nonconvex MINLP. However, it can also be re-written as a convex MINLP in a lifted space by using its epigraph formulation. To this end, we model each term in the objective function using a separate variable and bound it in a newly introduced constraint. The resulting problem then reads min x,c,η  p∈P j∈[k] ηpj (2a) s.t.η pj ≥p−cj2−Mp(1−xpj), p∈P,j∈[k],(2b)  j∈[k] xpj =1,p∈P,(2c) xpj ∈{0,1},p∈P,j∈[k],(2d) cj∈B,j∈[k],(2e) ηpj ≥0,p∈P,j∈[k],(2f) where Mpare sufficiently large numbers. The objective function is linear now and we obtain the additional quadratic and convex constraints in (2b). 123 Journal of Global Optimization (2023) 87:133–189 137 For every p∈P,Mpcan be chosen to be the maximum distance of pto any other point ˜p∈P. An overestimation can be easily computed via Mp=M=(u1−1)2+···+(ud−d)2, where iand uiare the componentwise bounds of the bounding box given above. Moreover, for a given cluster assignment x, an optimal choice for the cluster centroids is immediate, an observation that we will exploit frequently. Observation 2.1 For a given assignment of x-variables adhering to (1b)or (2c), respectively, the optimal choice for cj,j∈[k],is p∈Ppxpj p∈Pxpj , i.e., the barycenter of all points assigned to cluster j. 3 Cutting planes Without doubt, cutting planes are among the most powerful techniques to enhance the solution process for mixed-integer problems. Modern MI(N)LP solvers have many general-purpose cutting planes built-in. However, it is very often beneficial to derive problem-specific cutting planes. This is particularly important for the MSSC problem, since it is well-known that one of the most challenging issues for developing an efficient branch-and-bound algorithm for the MSSC problem is the computation of good lower bounds in a reasonable amount of time. In this section, we state two tailored cutting planes. The first one is applicable to Model (1) and (2) whereas the second one is only applicable to Model (2). 3.1 Cardinality cuts We first briefly discuss cardinality cuts, which are already mentioned in Aloise et al. [5]as well as Sherali and Desai [55]. Consider an optimal solution of Model (1). In this optimal solution, there cannot be any empty cluster, because otherwise the corresponding objective value can be decreased by assigning a point that is not a centroid to that empty cluster. On the other extreme, a cluster contains |P|−k+1 data points if every other cluster consists of only a single data point. Thus, to tighten the formulation (1), the following cardinality cuts can be added to Model (1): 1≤ p∈P xpj ≤|P|−k+1,j∈[k]. Obviously, the same cuts are also valid for Model (2). Moreover, note that the upper bound is implied by the model’s constraints and the lower bound of the previous inequalities as p∈Pk j=1xpj =|P|implies for a fixed j∈[k]that p∈Pxpj =|P|− p∈Pj∈[k]\{ j}xpj≤|P|−(k−1). The idea of cardinality cuts can also be localized, i.e., the cardinality bounds can be adapted to take local variable bounds at a node of the branch-and-bound tree into account. To this end, we introduce, for each j∈[k]the integral variable κjwith range {k−1,...,|P|−1} and link it with the x-variables via the linear constraint κj+p∈Pxpj =|P|,j∈[k]. 123 138 Journal of Global Optimization (2023) 87:133–189 That is, κjdescribes the number of data points that are not assigned to cluster j. If a lower bound κjand an upper bound ¯κjon κjis given, this equation implies the inequalities 1≤|P|−¯κj≤ p∈P xpj ≤|P|−κj≤|P|−k+1,j∈[k]. That is, they describe localized versions of cardinality cuts that get stronger if x-variables get fixed. Another side effect of the auxiliary variables κjis that a solver might decide to branch on these variables. In doing so, it imposes bounds on the size of cluster j∈[k]. 3.2 Outer approximation cuts We now focus on Model (2). The only nonlinear constraints (2b) in this problem are convex. Thus, their first-order Taylor approximations are global underestimators at any point (¯η, ¯c,¯x) and thus provide valid inequalities that are linear in (η, c,x): d  i=1 2¯cj icj i−2picj i+(pi)2−(¯cj i)2−ηpj −Mp(1−xpj)≤0,p∈P,j∈[k].(3) This allows to solve Model (2) in an outer approximation or LP/NLP-based branchand-bound fashion; see Duran and Grossmann [17], Fletcher and Leyffer [20]orQuesada and Grossmann [50], respectively. We start by relaxing the constraint set (2b). Next, we assume that (¯η, ¯c,¯x)is a solution of this relaxation, i.e., it particularly fulfills the binary conditions (2d). If the relaxation’s solution is feasible for the nonlinear constraints (2b), it is also a solution for Model (2). If not, we can compute a feasible point (ˆη, ˆc,¯x)of Model (2). In the original outer approximation method, this is done by fixing the binary variables ¯xin Model (2) and solving the resulting convex NLP subproblem; see Duran and Grossmann [17]. The benefit in our specific application is that solving the subproblem boils down to a simple computation of the barycenters ˆc, see Observation 2.1, followed by an evaluation of the distances ˆηaccording to Constraints (2b)and(2f). From the theory of outer approximation, it is well-known that when adding the inequalities (3) at the solution (ˆη, ˆc,¯x)of the subproblem, it holds  p∈P j∈[k] ηpj ≥ p∈P j∈[k] ˆηpj for all feasible points (η, c,¯x)of the updated relaxation. In other words, adding the outerapproximation cuts bounds the optimal objective value of the relaxation with fixed binaries x=¯xfrom below by p∈Pj∈[k]ˆηpj. Consequently, the updated relaxation yields a solution with a new, previously unseen, cluster assignment xor the optimality gap is closed. Thus, iterating this process terminates after a finite number of steps; see Duran and Grossmann [17] or Duran and Grossmann [20] for more details. We note that the number of inequalities (3) does not depend on the dimension d, which might be beneficial for problems with higher dimensions. Instead of implementing an LP/NLP-based branch-and-bound from scratch, we can use solvers such as SCIP to solve Model (2). In this setting, we can separate and add cuts (3)to tighten the LP relaxations. Since for larger |P|and k, adding all inequalities (3) might be impracticable, we may also add only a certain amount of cuts. In particular, in our implementation in SCIP, we add only 10 cuts per separation round. 123 Journal of Global Optimization (2023) 87:133–189 139 4 Propagation Suppose we are at a node of the branch-and-bound tree. Due to branching decisions and further reductions, some variables might have been fixed or their bounds have been tightened in comparison to the original problem formulation. The aim of propagation is to find further variable fixings or bound tightenings that are valid at the current node. That is, one tries to apply further reductions based on local variable bound information. According to Observation 2.1, every assignment of x-variables that satisfies (1b)or(2c) can be extended to a feasible solution of (1)or(2), respectively. Thus, it is crucial to derive propagation mechanisms that exclude assignments of x-variables that cannot be optimal. Moreover, we develop algorithms to strengthen bounds of the c-variables and the objective variables. Before we discuss our propagation algorithms, we fix the following notation and terminology. For every j∈[k],wedenotebyPj⊆Pthe set of all data points p∈Pwhose corresponding variable xpj has been fixed to 1 at the current node of the branch-and-bound tree. That is, we have already decided to assign pto cluster j. Moreover, we denote by P j⊆P all data points psuch that xpj has not been fixed to 0 yet, i.e., pis already or can still be assigned to cluster j. Note that Pj⊆P j. For a continuous variable z, i.e., for the cand η-variables,wedenotebyzand ¯zthe lower and upper bound on zat the current node, respectively. 4.1 Barycenter propagation Given a non-empty set of data points Q⊆Pdefining a cluster, the optimal choice for its centroid is the barycenter C(Q):= 1 |Q| p∈Q p of all data points in Q. The respective sum of all squared distances thus is D(Q)= p∈Qp−C(Q)2. The idea of the barycenter propagation is to use this observation to find lower bounds on the objective and to strengthen the bounds for the c-variables. 4.1.1 Bound tightening for the objective function values To find a lower bound on the objective in Model (1), note that for sets Q⊆Q⊆P, we have D(Q)≤D(Q). Consequently, a lower bound on the objective is given by j∈[k]D(Pj). The barycenter propagator uses this value to possibly tighten the lower bound on the objective at the current node of the branch-and-bound tree. Computing this lower bound for all clusters can be done in O(kd|P|)time and it has also been used by Brusco [7], see also Guns et al. [30], in a repetitive branch-and-bound framework. For Model (2), no immediate lower bound on the objective can be enforced because the objective is decoupled via the η-variables. Nevertheless, for each (p,j)∈P×[k],the following steps can be done. We can prune a node of the branch-and-bound tree if p/∈P j and ηpj >0, because an optimal solution has ηpj =0 as data points not assigned to a cluster do not contribute to D(Pj). Otherwise, if p/∈P jand ηpj =0, we can fix ηpj to 0. The first step is thus a pruning operation based on sub-optimal bounds in the subproblem, whereas the second step is a bound tightening operation. 123 140 Journal of Global Optimization (2023) 87:133–189 4.1.2 Bound tightening for the centroids Recall that [k]={1,...,k}. Besides strengthening bounds on the objective, barycenter information can also be used to tighten bounds on centroid variables cj iwith (i,j)∈[d]×[k]. Suppose P j\Pj={p1,...,ps}such that p1 i≤p2 i≤···≤ ps i. For each r∈[s]0:=[s]∪{0}, we compute γj,r=C(Pj∪{p1,...,pr}), i.e., the barycenter of the data points contained in Pjand the data points with the rsmallest ith coordinates that are not contained in Pj.As we show next, the ith coordinates of these barycenters can be used to compute a lower bound on cj i. Lemma 4.1 A valid lower bound on c j iis given by minr∈[s]0γj,r i. Proof Let Q⊆P j\Pjand assume |Q|=r. Then, C(Pj∪Q)i≥C(Pj∪{p1,...,pr})i, because the points p1,...,prare points with the rsmallest ith coordinates. Consequently, to find a lower bound on the centroids, it is sufficient to consider C(Pj∪{p1,...,pr})for each r∈[s]0. Analogously, an upper bound is given by maxr∈[s]0C(Pj∪{ps−r,...,ps})i. Since computing an iterative sequence of barycenters can be done using the formula γj,r+1=(|Pj|+r)γ j,r+pr+1 |Pj|+r+1, we can compute the minimum and maximum values for all coordinates and clusters in O(kd|P|)time. 4.2 Convexity and cone propagation Based on optimality arguments, we can also derive rules to assign data points p∈P j\Pjto cluster j∈[k]. The key idea of the convexity propagator is the following simple observation. Lemma 4.2 There exists an optimal solution of MSSC with clusters P1,...,Pksuch that, for each j ∈[k], we have conv(Pj)∩P=Pj. Proof Given an optimal allocation of the kcentroids, the Voronoi cells Cj={x∈Rd:x−cj≤x−cj,j∈[k]} for j∈[k]cover the entire Rdand only intersect at their boundaries. Since Voronoi cells are full-dimensional polyhedra, we can use the following mechanism to prove the assertion. We start with cluster 1 and observe that P1⊆C1in any optimal solution. If there exist p∈P\P1 that are contained in C1, they are necessarily contained in the boundary of C1. Hence, if we change the assignment of these points to P1, this does not change the objective of MSSC. The assertion thus holds for P1, and we can use the same arguments iteratively to conclude the proof.  As a consequence, the convexity propagator computes conv(Pj)for each j∈[k].Ifthere exists p∈P∩conv(Pj)it performs the following steps: If p/∈P jholds, then we can prune the current node of the branch-and-bound tree, because the local variable bounds cannot lead to an optimal solution adhering to Lemma 4.2. Otherwise, xpj can be fixed to 1. Besides pruning nodes and fixing variables to 1, Lemma 4.2 has another consequence that allows us to fix some variables to 0, which is illustrated in Fig. 1. 123 Journal of Global Optimization (2023) 87:133–189 147 To conduct the experiments, we use different test sets from the literature, which contain both real-world as well as synthetic instances. The test sets and the general computational setuparedescribedinSects.8.1 and 8.2, respectively. Then, in Sect. 8.3,westartthediscussion of the numerical results for the case k=2. We evaluate the benefits of each particular technique and indicate which setting performs best. Next, in Sect.8.4, we repeat the discussion but for the case k=3. Finally, in Sect.8.5, we present results on a larger test set in order to draw solid and comprehensive conclusions about the performance of the novel techniques. 8.1 Test sets We evaluate the impact of the presented algorithmic ideas for solving the MSSC problem using both synthetic and real-world test sets. To be able to draw conclusions on a reliable basis, we have collected all publicly available instances that have been used in the related literature for solving the MSSC problem to global optimality. Thus, to the best of our knowledge, our results are based on the largest publicly available test set for the MSSC problem consisting of realistic instances. Specifically, we use the instances that have been used in Aloise and Hansen [4], Sherali and Desai [55], as well as in Aloise et al. [5]. Since these instances come from different sources, we provide the source for every instance in Table 1. The synthetic test set has been proposed in Fränti and Sieranoja [22]. The authors show that these synthetic instances cover a wide range of classic MSSC instances. In particular, the test set contains instances with different degrees of overlap, density, and sparsity of data points. Note that some of the synthetic instances contain data points with very large coordinate values. In preliminary experiments, we have observed that this leads to very large big-M values in Model (2), which in turn causes numerical instabilities. To avoid numerical issues, we therefore re-scale these instances as follows. First, for each coordinate, we shift the data points such that their coordinate-wise minimum and maximum value is the same to have a “symmetric” distribution. Then, we re-scale the data points if they do not fit into [−103,103]d. More precisely, for each dimension i∈[d], we compute the maximum and minimum coordinate value obtaining ¯viand vi, respectively. Then, we take ui=0.5(¯vi+vi)and shift each data point pobtaining ˆpi=pi−uifor all i∈[d]. If we do this for all data points, they get centered around the origin. Now, if wi=vi−ui<−103or ¯wi=¯vi−ui>103holds, we re-scale the data. The desired new bounds then are zi=−103and ¯zi=103. Thus, the re-scaled data point ˜pis ˜pi=ˆpi−wi ¯wi−wi ·(¯zi−zi)+zi,i∈[d]. The corresponding instances that needed to be re-scaled are s1,s2,s3,s4,andunbalance. 8.2 Computational setup To conduct our experiments, we use SCIP 7.0.3 as a branch-and-bound framework. All LP relaxations are solved using CPLEX 12.8. Our novel techniques discussed in Sects.3–6are implemented as SCIP plugins written in C/C++ and our code is publicly available at GitHub2 (git hash 19003a37). To handle symmetries, we use the orbitope constraint handler plugin of SCIP, which implements orbitopal fixing and the symmetry handling inequalities as mentioned in Sect.7. To compute convex hulls and cones in the convexity propagator proposed in 2https://github.com/christopherhojny/globally-solving-MSSC. 123 148 Journal of Global Optimization (2023) 87:133–189 Table 1 Information about the test sets ID Instance Reference nd 1Fisher150iris Dua and Graff [16]andFisher[19] 150 4 2German22 Späth [57]222 3German59 Späth [57]592 4body-measurements Heinz et al. [33] 507 5 5cities-coord-202 Grötschel [29] 202 2 6cities-coord-666 Grötschel [29] 666 2 7concrete-compressive Dua and Graff [16] 1030 8 8glass-identification Dua and Graff [16] 214 9 9image-segmentation Dua and Graff [16] 2310 19 10 padberg-rinaldi-hole-dri Padberg and Rinaldi [42] 2392 2 11 reinelt-hole-drilling Reinelt [51] 1060 2 12 ruspini Ruspini [52]752 13 telugu-indian-vowel Pal and Majumder [43] 871 3 14 a1 Fränti and Sieranoja [22] 3000 2 15 a2 Fränti and Sieranoja [22] 5250 2 16 a3 Fränti and Sieranoja [22] 7500 2 17 dim Fränti and Sieranoja [22] 1024 32 18 g2-2-30 Fränti and Sieranoja [22] 2048 2 19 g2-2-50 Fränti and Sieranoja [22] 2048 2 20 g2-2-70 Fränti and Sieranoja [22] 2048 2 21 s1 Fränti and Sieranoja [22] 5000 2 22 s2 Fränti and Sieranoja [22] 5000 2 23 s3 Fränti and Sieranoja [22] 5000 2 24 s4 Fränti and Sieranoja [22] 5000 2 25 unbalance Fränti and Sieranoja [22] 6500 2 The first part corresponds to real-world test sets whose instances come from different sources. The second part corresponds to the synthetic test set Sect.4.2,weusetheQhull3C++ interface proposed by Barber et al. [6]. We have also conducted experiments using the CDD library [24] for computing convex hulls and cones, but due to numerical instabilities therein, we decided to use Qhull. Moreover, since we observed that many of SCIP’s internal heuristics require a lot of running time without generating a feasible solution, we disabled these heuristics. A list of disabled heuristics can be found in “Appendix A”. All computations were performed on a computer with two Intel Xeon CPU E5-2699 v4 at 2.20 GHz (2 ×44 threads) and 756 GB RAM. The time limit of all computations is 1h per instance. In the following, we discuss the impact of our techniques on solving the MSSC problem for k=2andk=3 clusters. We only report on aggregated results in the discussion and refer the reader to “Appendix B” for results per instance. The tables that we present show for both the quadratic and epigraph formulation the mean number of nodes in the branch-and-bound tree (column #nodes), the mean running time per instance in seconds (time), and the number of solved instances (#solved). Instances that cannot be solved within the time limit contribute 3http://www.qhull.org. 123 Journal of Global Optimization (2023) 87:133–189 149 Table 2 Comparison of mean number of nodes, mean running time per instance (in seconds), and number of solved instances using different heuristics for 2 clusters Setting Quadratic model Epigraph model round impr init #nodes time #solved #nodes time #solved 0 0 0 14,218.4 3406.45 1 23,275.1 2959.08 1 1 1 1 17,106.6 3273.37 1 19,884.9 2945.94 1 3600s to the mean time value. Moreover, we report on the used setting, where each of the following subsections describes how the settings are encoded in the tables. All mean numbers of measurements t1,...,tnare provided as shifted geometric means n i=1(ti+s)1 /n−sto reduce the impact of outliers. For time we use a shift of s=10 and for nodes a shift of s=100. 8.3 Discussion of the numerical results for 2 clusters We start with the discussion of the numerical results for the case when there are 2 clusters. First, we apply plain SCIP to all the 25 instances presented in Table 1. Afterward, we gradually enable our techniques in SCIP and evaluate the benefits of each particular technique as well as the benefits of different combinations of techniques. To allow for a concise encoding, we abbreviate the different techniques as described below. Whether a technique is enabled (resp. disabled) is encoded by 1 (resp. 0) in the corresponding tables. 8.3.1 Primal heuristics We start by evaluating the impact of primal heuristics. A summary of the obtained results is presented in Table 2, where “round.”, “impr.”, and “init”, serve as abbreviations for the rounding, improvement, and root-node heuristic, respectively. Recall that the improvement heuristic is only active for k>2, i.e., it has no effect in the experiments discussed next. The first row shows the results obtained by plain SCIP. It can be directly seen that the MSSC problem is extremely hard to solve. Note that SCIP is able to solve only 1 instance to global optimality, regardless of which model is used. Enabling all primal heuristics still does not allow to solve more instances. However, we can see that the mean running time decreases in both models, where the impact is larger for the quadratic model. That is, the single instance that can be solved is solved in approximately 3.9% faster using heuristics. Let us stress that, based on preliminary experiments, the main difficulty of solving the MSSC problem to global optimality is to obtain good dual bounds in a reasonable amount of time. For this reason, the impact of heuristics on the solving process is expected to be minor in comparison to the impact of techniques that improve the dual bound. However, since we needed to disable many of SCIP’s internal heuristics as described above, we enable all our heuristics in the following experiments as their running time is low and they produce good solutions. 8.3.2 Propagators The results of our experiments regarding propagators are summarized in Table 3,where “bary.”, “conv.”, “cone”, and “dist.” abbreviate the barycenter, convexity, cone, and distance 123 150 Journal of Global Optimization (2023) 87:133–189 Table 3 Comparison of mean number of nodes, mean running time per instance (in seconds), and number of solved instances using different propagators for 2 clusters Setting Quadratic model Epigraph model bary conv cone dist #nodes time #solved #nodes time #solved 0 0 0 0 17,106.6 3273.37 1 19,884.9 2945.94 1 1 0 0 0 91,401.7 2853.27 1 18,868.4 2863.78 2 0 1 0 0 3607.3 1735.34 4 10,478.6 1902.80 4 0 1 1 0 3029.8 1565.24 5 7448.3 1529.15 6 1 1 1 0 5996.1 835.39 8 6784.9 1527.04 5 propagator, respectively. However, we do not include the results using the distance propagator here, since in preliminary numerical experiments we observed that this propagator is not able to derive many reductions if used alone. In later experiments, we will enable it again to investigate whether it is able to improve the solution process if also other components are enabled. The first row of results in Table 3corresponds to the setting where only primal heuristics are enabled. It can be directly seen that as more propagators are enabled, more instances are solved. Without propagators only 1 instance is solved to global optimality. Using all our propagators, we are able to solve 8 instances with the quadratic model. Thus, the geometric ideas incorporated into the propagators are an important component to solve the MSSC problem effectively. In particular, plain SCIP is not able to make use of the simple geometric observations on its own. In the following, we discuss the benefits of each particular propagator in more detail. 8.3.3 Barycenter propagator By using the barycenter propagator and the quadratic model, much more nodes can be processed if compared with the previous setting and, more importantly, in significantly less time. The reason for this is that the barycenter propagator is able to perform many reductions, which in turn simplifies the LP relaxations. As a consequence, the dual bounds obtained with the quadratic model drastically improve by using the barycenter propagator. This can be clearly seen in Fig.2, where we plot the instances vs. the corresponding gap between the primal and dual bounds. This already demonstrates the great benefit that the barycenter propagator adds to the solution process. As discussed in Sect.4.1.1, the barycenter propagator is less powerful for the epigraph model as it is for the quadratic model, which is also reflected in the results. Nevertheless, it allows 1 more instance to be solved to global optimality and it slightly reduces the running times. Looking at Fig.2again, we also see that for many instances the gaps improve. 8.3.4 Convexity+Cone propagator The convexity propagator is based on geometric ideas and is extremely powerful. Using only this propagator alone and heuristics, we can already solve 3 more instances if the quadratic model is used, and 2 more instances if the epigraph model is used. Without the convexity propagator, these instances cannot be solved. This technique drastically helps in the solution 123 Journal of Global Optimization (2023) 87:133–189 151 Fig. 2 Instance ID vs. gap (in percentage and log-scale) for 2 clusters. Since the y-axis is in log-scale, if a particular instance is solved to global optimality by a particular method, then the gap is zero and hence it does not appear in the plot. Whereas, if the gap is equal or larger than 104, then it assumes the gap limit of 104in the plot process of the MSSC problem. Besides allowing more instances to be solved, it also requires half of the time that was needed before. Moreover, the number of nodes that need to be processed to solve the instances also reduces significantly. Using the cone propagation in combination with the convexity propagator, this effect is even more pronounced. It allows 1 more instance to be solved if the quadratic model is used, and 2 more instances if the epigraph model is used and results in much lower mean running times. 8.3.5 Barycenter+Convexity+Cone propagators Although the barycenter and convexity-cone propagators alone already significantly improved SCIP’s performance, their combination allows to solve three further instances in the quadratic model. This results in a significant reduction of running time by approximately 46%. Interestingly, the mean number of nodes in the combined setting is roughly twice as large as if just the convexity-cone propagator is used. This again shows that the reductions found by the propagators simplify the structure of relaxations drastically, e.g., because fixed x-variables remove non-convex expressions from the quadratic model. These reductions allow SCIP to process more nodes, which in turn allows to solve more instances. In the epigraph model, the combination of the three propagators does not qualitatively change the results. Although 1 less instance can be solved, for many instances the gaps improved; see Fig. 3. We conclude that, for both the quadratic and epigraph model, our propagation algorithms are an important component to solve the MSSC problem to global optimality. In particular, using combinations of these propagators creates synergies that allow to solve more instances in comparison with just using a single propagator, where the effect is more prominent for the quadratic model. 8.3.6 Cutting planes Next, we evaluate the impact of cutting planes on SCIP. As before, we also enable all heuristics and, due to the positive effect of propagators, also the convexity-cone and barycenter propagator. Preliminary numerical results showed that by localizing the cardinality cuts, no 123 152 Journal of Global Optimization (2023) 87:133–189 Fig. 3 Comparison of gaps using different propagators for 2 clusters Table 4 Comparison of mean number of nodes, mean running time per instance (in seconds), and number of solved instances using or not the OA cuts for 2 clusters Setting Epigraph model #nodes time #solved w/o OA cuts 6784.9 1527.04 5 w/ OA cuts 1526.8 1552.97 5 positive impact on the solution process can be achieved in general. Therefore, we focus only on the outer-approximation (OA) cuts. Since these cuts are only applicable for the epigraph model, we concentrate only on the epigraph model in the following discussion. Table 4 summarizes our results. At first glance, it seems that OA cuts only have a minor impact on SCIP’s performance as the number of solved instances does not change. Comparing the gaps with and without OA cuts, however, reveals a clear impact; see Fig.4. For 8 instances, we observe a change in the gap if cuts are enabled. In three cases, the gaps slightly degrade when OA cuts are enabled. For the remaining five instances, however, OA cuts either reduce or drastically reduce the gap. Thus, although no clear trend is visible, we may conclude that OA cuts are helpful when solving MSSC problems. The effect of OA cuts is less pronounced compared to the effect of propagators, which might be explained by the fact that OA cuts do not exploit the specific problem structure of MSSC. In contrast to this, our novel propagator techniques are tailored to the MSSC problem and thus allow stronger reductions. 8.3.7 Distance propagator As reported above, the distance propagator alone is not able to significantly improve SCIP’s performance. For this reason, we test its effect if also further components are enabled. From Table 5, we can see that using the distance propagator together with our other techniques has a slightly positive effect. Therefore, it is also enabled in the experiments discussed next. 123 Journal of Global Optimization (2023) 87:133–189 153 Fig. 4 Comparison of gaps using propagators with the OA cuts enabled or not for 2 clusters Table 5 Comparison of mean number of nodes, mean running time per instance (in seconds), and number of solved instances using or not distance propagator for 2 clusters Setting Quadratic model Epigraph model #nodes time #solved #nodes time #solved w/o dist. propagator 5996.1 835.39 8 1526.8 1552.97 5 w/ dist. propagator 5976.4 813.00 8 1809.9 1497.49 6 8.3.8 Branching rules The last components to be tested are branching rules. We have implemented all branching rules described in Sect.5. Preliminary numerical experiments, however, revealed that only the entropy and distance branching rules may be beneficial for some instances. In contrast, the centrality and pairs-in-the-middle-of-pairs harm the solution process leading to a less well-performing code. Therefore, we focus only on the branching rules that have a positive impact on some instances in the following discussion. We present the summary results in Table 6.By“standard”werefertoSCIP’s default branching rule. Our experiments show that no branching rule dominates the others. On the one hand, by using the distance branching rule and the quadratic model, 1 additional instance can be solved. On the other hand, the number of nodes and the time required increases. If the epigraph model is used, then the entropy branching rule is performing best: The number of solved instances remains the same but the running times are slightly lower. The overall impact of branching rules, however, seems to heavily depend on the underlying instance to solve, which does not allow us to provide a clear winner. 123 154 Journal of Global Optimization (2023) 87:133–189 Table 6 Comparison of mean number of nodes, mean running time per instance (in seconds), and number of solved instances using different branching rules for 2 clusters Setting Quadratic model Epigraph model #nodes time #solved #nodes time #solved Standard 5976.4 813.00 8 1809.9 1497.49 6 Entropy 6613.7 863.79 8 2541.1 1446.77 6 Distance 6942.0 845.73 9 1878.8 1488.07 6 Fig. 5 Running times and gaps comparison between plain SCIP and SCIP enabled with the best setting for 2 clusters 8.3.9 Best setting To conclude the discussion of the numerical results for k=2, we show a comparison of plain SCIP with the best combination of the techniques proposed in this paper. The latter comprises primal heuristics, propagators, OA cuts (for the epigraph model), and the standard branching rules of SCIP, since our branching rules and standard branching rules are performing equally good on average. This comparison in shown in Fig. 5. Regarding the performance of plain SCIP, we emphasize that the dual bounds found by SCIP in the quadratic model are very weak which leads to very large gaps. In contrast to this, if the epigraph formulation is used, better dual bounds can be obtained by using plain SCIP. Despite their simplicity, the plots show that our novel geometric ideas drastically improve 123 Journal of Global Optimization (2023) 87:133–189 155 Table 7 Comparison of mean number of nodes, mean running time per instance (in seconds), and number of solved instances using different heuristics for 3 clusters Setting Quadratic model Epigraph model round impr init #nodes time #solved #nodes time #solved 0 0 0 7365.1 3600.00 0 5756.0 3094.27 1 1 1 1 8400.1 3600.00 0 6998.6 3113.39 1 on the performance of SCIP, thus adding powerful methods to the toolbox for solving the MSSC problem to global optimality if k=2. These methods work particularly well for instances with 2-dimensional data and the number of data points not exceeding 2048 as almost all such instances from our test set, see Table 1, can be solved to global optimality within the time limit by the quadratic model. The only exception is instance 11, which terminates after one hour with a gap of 19.02%. For more detailed results of the best setting we refer the reader to Table 21 in the “Appendix B”. 8.4 Discussion of the numerical results for 3 clusters We now turn our attention to the experiments for k=3. The MSSC problem is much harder to solve to global optimality in this setting. We proceed as in the last section. 8.4.1 Primal heuristics The summary results of plain SCIP and SCIP enabled with our primal heuristics are presented in Table 7. By using plain SCIP and the epigraph model, we can solve only 1 instance to global optimality. Enabling heuristics does not change the number of solved instances. This is in line with the observations made above: the main difficulty in solving the MSSC problem to global optimality is to provide tight dual bounds, which are not provided by primal heuristics. However, we observe that by enabling the heuristics in the epigraph model, the primal-dual gap improves for many instances substantially; see Fig.6. For this reason, we enable heuristics in the following experiments. 8.4.2 Propagators Next, we evaluate the impact of our propagation techniques for k=3. In Table 8, we show the summarized results obtained by activating our primal heuristics and the propagators. Taking a general look at the results and comparing the first row (SCIP + heuristics) with the last row, we can see that 2 more instances can be solved to global optimality, using either the quadratic or the epigraph model. Thus, although the MSSC problem for k=3ismuchhardertosolve than for k=2, the propagation techniques are still helpful; see also Fig. 7where we compare the gaps obtained by using propagators. Therefore, in what follows, we discuss the benefits of the separate propagators in turn. 123 156 Journal of Global Optimization (2023) 87:133–189 Fig. 6 Comparison of gaps using or not the heuristics for 3 clusters Table 8 Comparison of mean number of nodes, mean running time per instance (in seconds), and number of solved instances using different propagators for 3 clusters Setting Quadratic model Epigraph model bary conv cone dist #nodes time #solved #nodes time #solved 0 0 0 0 7365.1 3600.00 0 5756.0 3094.27 1 1 0 0 0 20,902.8 2983.99 1 7409.5 2977.89 1 0 1 0 0 8066.6 3250.65 1 9169.4 2939.98 1 0 1 1 0 11,815.9 3385.49 1 7196.0 2743.45 3 1 1 1 0 26,530.7 2856.07 2 6057.0 2760.58 3 Fig. 7 Comparison of gaps using different propagators for 3 clusters 123 Journal of Global Optimization (2023) 87:133–189 163 Table 12 continued Setting Quadratic model Epigraph model bary conv cone dist #nodes time #solved #nodes time #solved Sample (600) 0 0 0 0 9475.0 3600.00 0 29,099.4 3600.00 0 1 0 0 0 181,258.9 3600.00 0 21,186.1 3600.00 0 0 1 0 0 28,615.5 3600.00 0 28,174.6 3600.00 0 0 1 1 0 21,447.4 1941.54 4 7426.4 913.55 7 1 1 1 0 5837.4 61.37 6 8077.5 948.55 7 1 1 1 1 5863.2 61.20 6 8077.5 974.16 7 Sample (700) 0 0 0 0 7734.4 3600.00 0 26 465.3 3600.00 0 1 0 0 0 152,656.8 3600.00 0 19,163.8 3600.00 0 0 1 0 0 23,220.9 3600.00 0 21,774.4 3600.00 0 0 1 1 0 30,603.9 1780.14 6 7814.5 1131.49 7 1 1 1 0 6372.9 138.53 5 8800.0 1257.95 7 1 1 1 1 6325.5 140.55 5 8800.0 1274.35 7 Sample (800) 0 0 0 0 6139.6 3600.00 0 22,311.9 3600.00 0 1 0 0 0 130,506.2 3600.00 0 17,108.1 3600.00 0 0 1 0 0 18,217.5 3600.00 0 17,203.4 3600.00 0 0 1 1 0 26,364.9 2506.24 4 9310.2 1518.02 7 1 1 1 0 7861.5 80.38 6 9628.0 1570.32 7 1 1 1 1 7840.7 79.63 6 9628.0 1546.33 7 Sample (900) 0 0 0 0 5808.0 3600.00 0 20,471.8 3600.00 0 1 0 0 0 118,504.8 3600.00 0 15,074.9 3600.00 0 0 1 0 0 16,275.8 3600.00 0 13,287.0 3600.00 0 0 1 1 0 25,767.8 3174.30 3 9624.2 1885.05 7 1 1 1 0 5741.3 301.63 4 9732.6 1898.40 7 1 1 1 1 7102.9 177.60 5 9732.6 1862.99 7 Sample (1000) 0 0 0 0 4839.0 3600.00 0 17,887.6 3600.00 0 1 0 0 0 104,220.5 3600.00 0 14,923.6 3600.00 0 0 1 0 0 12,403.4 3600.00 0 11,910.6 3600.00 0 0 1 1 0 23,446.3 3372.88 2 10,227.0 2127.64 6 1 1 1 0 8852.9 109.65 6 10,752.2 2148.66 6 1 1 1 1 7873.6 182.62 5 10,840.9 2089.13 6 123 164 Journal of Global Optimization (2023) 87:133–189 Table 13 Comparison of mean number of nodes, mean running time per instance (in seconds), and number of solved instances using or not the OA cuts for the epigraph model Epigraph model Setting #nodes time #solved Sample (100) w/o OA cuts 2171.7 42.65 7 w/ OA cuts 2459.9 50.21 7 Sample (200) w/o OA cuts 3094.1 124.43 7 w/ OA cuts 3327.1 177.07 7 Sample (300) w/o OA cuts 4387.3 262.19 7 w/ OA cuts 5051.0 385.62 7 Sample (400) w/o OA cuts 5430.1 408.90 7 w/ OA cuts 6171.4 624.12 7 Sample (500) w/o OA cuts 6998.6 681.69 7 w/ OA cuts 7249.4 1004.41 7 Sample (600) w/o OA cuts 8077.5 974.16 7 w/ OA cuts 7792.4 1292.30 7 Sample (700) w/o OA cuts 8800.0 1274.35 7 w/ OA cuts 9289.0 1754.04 7 Sample (800) w/o OA cuts 9628.0 1546.33 7 w/ OA cuts 10,139.6 2224.56 7 Sample (900) w/o OA cuts 9732.6 1862.99 7 w/ OA cuts 10,387.2 2300.25 6 Sample (1000) w/o OA cuts 10,840.9 2089.13 6 w/ OA cuts 9838.3 2622.94 4 From this experiment, we also conclude that, in general, the quadratic model performs better regarding running times, whereas the epigraph model performs better regarding the dual bounds, which in turn leads to more instances being solved. 9 Conclusion Solving the MSSC problem to global optimality is a very challenging task that already has received considerable attention in the literature. Nevertheless, the problem is far from being 123 Journal of Global Optimization (2023) 87:133–189 165 Table 14 Comparison of mean number of nodes, mean running time per instance (in seconds), and number of solved instances using different branching rules Setting Quadratic model Epigraph model #nodes time #solved #nodes time #solved Sample (100) Standard 1715.9 1.91 7 2459.9 50.21 7 Distance 2220.6 2.50 7 1990.6 46.67 7 Entropy 2133.9 2.38 7 2043.8 41.30 7 Sample (200) Standard 3229.1 6.78 7 3327.1 177.07 7 Distance 3715.7 8.33 7 3626.5 218.33 7 Entropy 3652.8 7.94 7 3572.1 164.08 7 Sample (300) Standard 4632.5 12.25 7 5051.0 385.62 7 Distance 5134.7 14.29 7 5174.6 493.21 7 Entropy 4445.0 35.91 6 4959.5 354.50 7 Sample (400) Standard 5356.3 16.11 7 6171.4 624.12 7 Distance 5601.1 18.33 7 5934.1 792.31 7 Entropy 5169.0 44.84 6 6164.3 571.01 7 Sample (500) Standard 5956.7 21.09 7 7249.4 1004.41 7 Distance 5086.9 58.55 6 7653.5 1182.79 7 Entropy 6205.8 59.10 6 7124.9 902.84 7 Sample (600) Standard 5863.2 61.93 6 7792.4 1292.30 7 Distance 7754.0 34.30 7 9088.0 1763.90 7 Entropy 6690.4 27.86 7 7922.3 1083.07 7 Sample (700) Standard 6325.5 139.71 5 9289.0 1754.04 7 Distance 7288.0 82.41 6 10,056.9 2147.21 7 Entropy 7625.7 34.99 7 9250.8 1540.31 7 Sample (800) Standard 7840.7 78.82 6 10,139.6 2224.56 7 Distance 6844.6 163.07 5 10,594.9 2518.07 5 Entropy 7535.0 151.02 5 10,435.5 1969.20 6 Sample (900) Standard 7102.9 174.98 5 10 387.2 2300.25 6 Distance 6872.1 166.80 5 8988.8 2390.23 5 Entropy 6307.4 158.35 5 11,012.2 2308.25 5 123 166 Journal of Global Optimization (2023) 87:133–189 Table 14 continued Setting Quadratic model Epigraph model #nodes time #solved #nodes time #solved Sample (1000) Standard 7873.6 180.79 5 9838.3 2622.94 4 Distance 9924.7 62.60 7 9021.6 2705.16 4 Entropy 8833.3 200.34 5 11,400.1 2663.63 5 Fig. 12 Running times and gaps comparison between plain SCIP and SCIP enabled with the best setting for the sampled instances and for 2 clusters “practically solved”. In this paper, we propose different techniques (including propagation, cutting planes, branching rules, or primal heuristics) that can be incorporated in a branch-andbound framework for solving the problem. Our extensive numerical study shows that these novel techniques significantly help to improve the solution process. On the one hand, we can now solve instances that have not been solvable before. On the other hand, the optimality gaps for those instances that remain unsolvable are significantly reduced. Not surprisingly, there are still some ideas left for future research. Let us sketch two of them. First, we show that our techniques can be used to globally solve instances of moderate size. Thus, our methods could also be used in solution approaches for the MSSC problem that rely on reducing the dimension or the size of the originally given problem; see, e.g., Hua et al. [35]. Second, there further exist variants of the MSSC problem with additional side constraints as discussed in, e.g., Liberti and Manca [38]. Such side constraints allow 123 Journal of Global Optimization (2023) 87:133–189 167 for solution techniques that are feasibility-based, whereas all our techniques are optimalitybased. Hence, a combination of both could yield an overall branch-and-bound framework that is even more effective for side-constrained MSSC problems. Acknowledgements The second author thanks the DFG for their support within RTG 2126 “Algorithmic Optimization”. Moreover, the last author thanks the DFG for their support within projects A05 and B08 in CRCTRR154. Funding Open Access funding enabled and organized by Projekt DEAL. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. Appendix A. List of SCIP heuristics disabled in our experiments We disabled the following SCIP heuristics: alns, mpec, subnlp, rens, locks, objpscostdiving, distributiondiving, clique, nlpdiving, rins, linesearchdiving, conflictdiving, crossover, fracdiving, guideddiving, pscostdiving, randrounding, veclendiving, adaptivediving. Appendix B. Numerical results per instance Here we present the detailed numerical results obtained for each instance. The running time (in seconds) is encoded as “time”. The gap (in percentage) between the primal and dual bounds is encoded as “gap”. See Tables 15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,and33. 123 168 Journal of Global Optimization (2023) 87:133–189 Table 15 Results for 2 clusters using plain SCIP Instance Quadratic model Epigraph model #nodes time gap #nodes time gap Fisher150iris 909,913 3600.00 ∞18,956 3600.00 170.24 German22 207,026 900.31 0.00 1343 17.25 0.00 German59 203,438 3600.00 ∞112,482 3600.00 32.98 body-measurements 38,934 3600.00 ∞8862 3600.00 3338.16 cities-coord-202 84,088 3600.00 ∞30,276 3600.00 284.64 cities-coord-666 16,459 3600.00 ∞19,784 3600.00 2886.05 concrete-compressive 12,496 3600.00 ∞8505 3600.00 79149.66 glass-identification 240,834 3600.00 ∞12,176 3600.00 1920.57 image-segmentation 3026 3600.00 ∞2870 3600.00 8405.89 padberg-rinaldi-hole-dri 2654 3600.00 ∞114,594 3600.00 ∞ reinelt-hole-drilling 2265 3600.00 ∞235,206 3600.00 ∞ ruspini 874,673 3600.00 ∞137,981 3600.00 132.28 telugu-indian-vowel 846 3600.00 ∞13,681 3600.00 5082.69 a1 3024 3600.00 ∞88,111 3600.00 ∞ a2 5338 3600.00 ∞67,159 3600.00 ∞ a3 6852 3600.00 ∞41,941 3600.00 ∞ dim 3175 3600.00 ∞738 3600.00 ∞ g2-2-30 13,927 3600.00 ∞10,273 3600.00 22444.38 g2-2-50 11,804 3600.00 ∞10,837 3600.00 13433.08 g2-2-70 2376 3600.00 ∞8235 3600.00 19192.83 s1 5032 3600.00 ∞58,633 3600.00 ∞ s2 5050 3600.00 ∞47,121 3600.00 ∞ s3 5035 3600.00 ∞55,728 3600.00 ∞ s4 5029 3600.00 ∞53,218 3600.00 ∞ unbalance 6544 3600.00 ∞37,468 3600.00 ∞ 123 Journal of Global Optimization (2023) 87:133–189 169 Table 16 Results for 2 clusters using enabled heuristics Instance Quadratic model Epigraph model #nodes time gap #nodes time gap Fisher150iris 818,860 3600.00 ∞15,462 3600.00 219.06 German22 112,049 327.16 0.00 1141 14.39 0.00 German59 266,506 3600.00 ∞142,167 3600.00 32.20 body-measurements 100,130 3600.00 ∞7959 3600.00 2441.41 cities-coord-202 74,901 3600.00 ∞37,155 3600.00 234.43 cities-coord-666 13,015 3600.00 ∞16,950 3600.00 1980.38 concrete-compressive 26,517 3600.00 ∞4933 3600.00 17012.20 glass-identification 253,985 3600.00 ∞10,690 3600.00 905.82 image-segmentation 2666 3600.00 ∞1601 3600.00 28016.73 padberg-rinaldi-hole-dri 4703 3600.00 ∞101,082 3600.00 ∞ reinelt-hole-drilling 2081 3600.00 ∞189,387 3600.00 ∞ ruspini 736,505 3600.00 ∞142,065 3600.00 85.06 telugu-indian-vowel 852 3600.00 ∞10,886 3600.00 2774.03 a1 5913 3600.00 ∞76,717 3600.00 ∞ a2 5156 3600.00 ∞24,795 3600.00 ∞ a3 7083 3600.00 ∞38,235 3600.00 ∞ dim 3991 3600.00 ∞1257 3600.00 ∞ g2-2-30 13,412 3600.00 ∞9103 3600.00 5636.36 g2-2-50 10,475 3600.00 ∞9832 3600.00 8397.81 g2-2-70 11,617 3600.00 ∞11,084 3600.00 8115.49 s1 4840 3600.00 ∞29,981 3600.00 ∞ s2 9902 3600.00 ∞47,364 3600.00 ∞ s3 9786 3600.00 ∞44,516 3600.00 ∞ s4 4917 3600.00 ∞45,528 3600.00 ∞ unbalance 6075 3600.00 ∞28,532 3600.00 ∞ 123 170 Journal of Global Optimization (2023) 87:133–189 Table 17 Results for 2 clusters using enabled heuristics and barycenter propagator Instance Quadratic model Epigraph model #nodes time gap #nodes time gap Fisher150iris 1,291,451 3600.00 390.75 36,671 3600.00 133.79 German22 735 1.00 0.00 603 3.14 0.00 German59 298,474 3600.00 60.67 181,024 3302.02 0.00 body-measurements 175,149 3600.00 4020.08 10,707 3600.00 2023.87 cities-coord-202 209,526 3600.00 182.11 53,429 3600.00 135.06 cities-coord-666 6326 3600.00 2030.10 29,584 3600.00 1624.17 concrete-compressive 78,715 3600.00 5124.69 5448 3600.00 18755.18 glass-identification 378,677 3600.00 327.75 13,634 3600.00 449.85 image-segmentation 13,938 3600.00 41771.35 251 3600.00 ∞ padberg-rinaldi-hole-dri 215,038 3600.00 5111.74 115,222 3600.00 ∞ reinelt-hole-drilling 548,760 3600.00 1821.46 243,529 3600.00 ∞ ruspini 371,336 3600.00 172.24 202,409 3600.00 62.62 telugu-indian-vowel 482,986 3600.00 1188.41 11,452 3600.00 2802.01 a1 154,455 3600.00 9892.16 76,186 3600.00 ∞ a2 79,802 3600.00 14609.94 28,009 3600.00 ∞ a3 50,747 3600.00 28536.14 8957 3600.00 ∞ dim 20,279 3600.00 12806.54 938 3600.00 ∞ g2-2-30 68,620 3600.00 10686.30 9677 3600.00 8018.12 g2-2-50 69,035 3600.00 8460.72 11,622 3600.00 6694.39 g2-2-70 75,071 3600.00 6784.35 10,224 3600.00 6502.58 s1 78,784 3600.00 29879.42 30,042 3600.00 ∞ s2 82,041 3600.00 18242.68 34,579 3600.00 ∞ s3 87,351 3600.00 17978.65 32,344 3600.00 ∞ s4 87,945 3600.00 17583.74 35,208 3600.00 ∞ unbalance 70,452 3600.00 5680.62 22,020 3600.00 ∞ 123 Journal of Global Optimization (2023) 87:133–189 171 Table 18 Results for 2 clusters using enabled heuristics and convexity propagator Instance Quadratic model Epigraph model #nodes time gap #nodes time gap Fisher150iris 1,262,448 3600.00 ∞38,295 3600.00 130.60 German22 621 1.47 0.00 523 2.55 0.00 German59 2887 10.46 0.00 1979 17.74 0.00 body-measurements 232,593 3600.00 ∞12,806 3600.00 1701.41 cities-coord-202 150,793 293.82 0.00 37,339 1333.65 0.00 cities-coord-666 179,651 3600.00 12.50 17,055 3600.00 93.41 concrete-compressive 771 3600.00 ∞6482 3600.00 31862.98 glass-identification 210,831 3600.00 ∞12,077 3600.00 885.51 image-segmentation 1946 3600.00 ∞2394 3600.00 28016.73 padberg-rinaldi-hole-dri 111 3600.00 ∞33,732 3600.00 ∞ reinelt-hole-drilling 231 3600.00 ∞131,094 3600.00 ∞ ruspini 9953 20.64 0.00 1957 36.14 0.00 telugu-indian-vowel 171 3600.00 ∞8986 3600.00 354.14 a1 183 3600.00 ∞40,666 3600.00 ∞ a2 377 3600.00 ∞20,874 3600.00 ∞ a3 64 3600.00 ∞12,275 3600.00 ∞ dim 3278 3600.00 ∞1502 3600.00 ∞ g2-2-30 107,172 3600.00 ∞5668 3600.00 29.91 g2-2-50 196,144 3600.00 ∞3547 3600.00 34.49 g2-2-70 226,646 3600.00 ∞2775 3600.00 92.31 s1 84 3600.00 ∞22,269 3600.00 ∞ s2 182 3600.00 ∞18,286 3600.00 ∞ s3 89 3600.00 ∞23,013 3600.00 ∞ s4 98 3600.00 ∞22,590 3600.00 ∞ unbalance 98 3600.00 ∞26,895 3600.00 ∞ 123 172 Journal of Global Optimization (2023) 87:133–189 Table 19 Results for 2 clusters using enabled heuristics and convexity+cone propagator Instance Quadratic model Epigraph model #nodes time gap #nodes time gap Fisher150iris 1,254,550 3600.00 ∞38,867 3600.00 114.71 German22 539 1.40 0.00 269 1.28 0.00 German59 2139 7.71 0.00 627 7.45 0.00 body-measurements 247,712 3600.00 ∞10,235 3600.00 1460.22 cities-coord-202 49,373 142.42 0.00 4759 155.58 0.00 cities-coord-666 63,323 773.03 0.00 3644 525.23 0.00 concrete-compressive 962 3600.00 ∞6464 3600.00 31862.98 glass-identification 216,426 3600.00 ∞12,824 3600.00 877.58 image-segmentation 2175 3600.00 ∞2193 3600.00 28016.73 padberg-rinaldi-hole-dri 53 3600.00 ∞44,395 3600.00 ∞ reinelt-hole-drilling 26 3600.00 ∞161,496 3600.00 ∞ ruspini 6875 15.21 0.00 549 10.29 0.00 telugu-indian-vowel 361 3600.00 ∞8627 3600.00 171.65 a1 340 3600.00 ∞40,987 3600.00 ∞ a2 22 3600.00 ∞10,912 3600.00 ∞ a3 219 3600.00 ∞5419 3600.00 ∞ dim 3318 3600.00 ∞1655 3600.00 ∞ g2-2-30 172,536 3600.00 ∞6249 3460.79 0.00 g2-2-50 205,102 3600.00 ∞4509 3600.00 11.25 g2-2-70 190,472 3600.00 ∞3714 3600.00 44.85 s1 206 3600.00 ∞15,937 3600.00 ∞ s2 58 3600.00 ∞14,702 3600.00 ∞ s3 70 3600.00 ∞14,203 3600.00 ∞ s4 34 3600.00 ∞16,303 3600.00 ∞ unbalance 24 3600.00 ∞20,735 3600.00 ∞ 123 Journal of Global Optimization (2023) 87:133–189 179 Table 26 Results for 3 clusters using enabled heuristics and barycenter propagator Instance Quadratic model Epigraph model #nodes time gap #nodes time gap Fisher150iris 280,605 3600.00 1734.55 28,221 3600.00 880.33 German22 15,992 23.58 0.00 6493 21.91 0.00 German59 314,057 3600.00 241.47 268,192 3600.00 251.06 body-measurements 43,177 3600.00 11222.63 8046 3600.00 ∞ cities-coord-202 103,462 3600.00 943.44 60,570 3600.00 423.04 cities-coord-666 11,930 3600.00 5752.10 19,756 3600.00 4775.80 concrete-compressive 6277 3600.00 ∞2448 3600.00 ∞ glass-identification 118,483 3600.00 1069.56 14,447 3600.00 1443.69 image-segmentation 3574 3600.00 ∞322 3600.00 ∞ padberg-rinaldi-hole-dri 78,593 3600.00 ∞450 3600.00 ∞ reinelt-hole-drilling 255,437 3600.00 833137.03 443,947 3600.00 ∞ ruspini 292,533 3600.00 715.82 211,537 3600.00 237.44 telugu-indian-vowel 230,083 3600.00 143344.25 6129 3600.00 30813.12 a1 59,081 3600.00 ∞142 3600.00 ∞ a2 11,561 3600.00 ∞594 3600.00 ∞ a3 0 3600.00 ∞74 3600.00 ∞ dim 6135 3600.00 ∞323 3600.00 ∞ g2-2-30 11,274 3600.00 399488.15 4714 3600.00 788662.43 g2-2-50 13,451 3600.00 206676.83 5048 3600.00 449014.04 g2-2-70 5850 3600.00 ∞4677 3600.00 ∞ s1 16,290 3600.00 ∞23,627 3600.00 ∞ s2 15,558 3600.00 ∞76,626 3600.00 ∞ s3 19,054 3600.00 ∞1999 3600.00 ∞ s4 20,110 3600.00 ∞35,733 3600.00 ∞ unbalance 1 3600.00 ∞15,560 3600.00 ∞ 123 180 Journal of Global Optimization (2023) 87:133–189 Table 27 Results for 3 clusters using enabled heuristics and convexity propagator Instance Quadratic model Epigraph model #nodes time gap #nodes time gap Fisher150iris 442,269 3600.00 ∞36,077 3600.00 497.90 German22 122,884 273.42 0.00 4506 13.19 0.00 German59 694,633 3600.00 ∞371,493 3600.00 58.32 body-measurements 60,146 3600.00 ∞7474 3600.00 14488.38 cities-coord-202 636,337 3600.00 ∞102,010 3600.00 737.75 cities-coord-666 26,103 3600.00 ∞45,587 3600.00 6491.97 concrete-compressive 1 3600.00 ∞1 3600.00 ∞ glass-identification 51,612 3600.00 ∞11,479 3600.00 9409.18 image-segmentation 977 3600.00 ∞388 3600.00 ∞ padberg-rinaldi-hole-dri 3483 3600.00 ∞2883 3600.00 ∞ reinelt-hole-drilling 1372 3600.00 ∞415,922 3600.00 ∞ ruspini 753,684 3600.00 ∞323,701 3600.00 248.52 telugu-indian-vowel 410,343 3600.00 ∞8846 3600.00 20873.74 a1 368 3600.00 ∞708 3600.00 ∞ a2 2688 3600.00 ∞942 3600.00 ∞ a3 0 3600.00 ∞74 3600.00 ∞ dim 1343 3600.00 ∞429 3600.00 ∞ g2-2-30 8084 3600.00 ∞14,399 3600.00 ∞ g2-2-50 12,276 3600.00 ∞5698 3600.00 104726.70 g2-2-70 3745 3600.00 ∞3829 3600.00 96142.33 s1 3447 3600.00 ∞53,673 3600.00 ∞ s2 582 3600.00 ∞24,445 3600.00 ∞ s3 322 3600.00 ∞11,619 3600.00 ∞ s4 448 3600.00 ∞49,473 3600.00 ∞ unbalance 4324 3600.00 ∞16,934 3600.00 ∞ 123 Journal of Global Optimization (2023) 87:133–189 181 Table 28 Results for 3 clusters using enabled heuristics and convexity+cone propagator Instance Quadratic model Epigraph model #nodes time gap #nodes time gap Fisher150iris 229,614 3600.00 ∞34,610 3600.00 531.20 German22 234,561 770.50 0.00 6162 21.35 0.00 German59 501,525 3600.00 ∞145,215 1285.61 0.00 body-measurements 55,863 3600.00 ∞8049 3600.00 12590.36 cities-coord-202 233,440 3600.00 ∞91,917 3600.00 275.14 cities-coord-666 26,776 3600.00 ∞18,046 3600.00 2763.85 concrete-compressive 1 3600.00 ∞1 3600.00 ∞ glass-identification 72,842 3600.00 ∞11,375 3600.00 9515.46 image-segmentation 882 3600.00 ∞384 3600.00 ∞ padberg-rinaldi-hole-dri 8466 3600.00 ∞2883 3600.00 ∞ reinelt-hole-drilling 1022 3600.00 ∞350,191 3600.00 ∞ ruspini 710,135 3600.00 ∞123,603 1317.54 0.00 telugu-indian-vowel 90,336 3600.00 ∞7121 3600.00 37240.30 a1 0 3600.00 ∞708 3600.00 ∞ a2 27,022 3600.00 ∞942 3600.00 ∞ a3 0 3600.00 ∞74 3600.00 ∞ dim 1579 3600.00 ∞395 3600.00 ∞ g2-2-30 23,034 3600.00 ∞4736 3600.00 440006.61 g2-2-50 21,069 3600.00 ∞4186 3600.00 7295513.20 g2-2-70 30,670 3600.00 ∞3783 3600.00 ∞ s1 10,659 3600.00 ∞49,155 3600.00 ∞ s2 1772 3600.00 ∞22,270 3600.00 ∞ s3 7235 3600.00 ∞2985 3600.00 ∞ s4 4709 3600.00 ∞49,473 3600.00 ∞ unbalance 3769 3600.00 ∞16,966 3600.00 ∞ 123 182 Journal of Global Optimization (2023) 87:133–189 Table 29 Results for 3 clusters using enabled heuristics and barycenter+convexity+cone propagator Instance Quadratic model Epigraph model #nodes time gap #nodes time gap Fisher150iris 348,209 3600.00 1890.73 29,323 3600.00 767.07 German22 7569 10.77 0.00 4576 18.24 0.00 German59 534,709 1949.13 0.00 143,355 1226.50 0.00 body-measurements 74,701 3600.00 15297.12 8419 3600.00 16458.64 cities-coord-202 218,931 3600.00 727.52 92,471 3600.00 246.47 cities-coord-666 26,602 3600.00 58814.95 18,797 3600.00 3456.63 concrete-compressive 1 3600.00 ∞1 3600.00 ∞ glass-identification 145,443 3600.00 1034.68 14,558 3600.00 1439.54 image-segmentation 4374 3600.00 ∞317 3600.00 ∞ padberg-rinaldi-hole-dri 150,154 3600.00 ∞450 3600.00 ∞ reinelt-hole-drilling 379,375 3600.00 ∞363,311 3600.00 ∞ ruspini 819,053 3600.00 871.99 145,815 1793.26 0.00 telugu-indian-vowel 281,326 3600.00 ∞9796 3600.00 637541.95 a1 111,630 3600.00 ∞142 3600.00 ∞ a2 30,873 3600.00 ∞594 3600.00 ∞ a3 0 3600.00 ∞74 3600.00 ∞ dim 7302 3600.00 ∞294 3600.00 ∞ g2-2-30 35,125 3600.00 99442.00 3825 3600.00 ∞ g2-2-50 34,063 3600.00 69060.41 3766 3600.00 960564.45 g2-2-70 12,876 3600.00 ∞5841 3600.00 293061.36 s1 5624 3600.00 ∞21,594 3600.00 ∞ s2 5496 3600.00 ∞65,969 3600.00 ∞ s3 1599 3600.00 ∞1365 3600.00 ∞ s4 65,193 3600.00 ∞35,733 3600.00 ∞ unbalance 5514 3600.00 ∞13,148 3600.00 ∞ 123 Journal of Global Optimization (2023) 87:133–189 183 Table 30 Results for 3 clusters using enabled heuristics, barycenter+convexity+cone propagators, and OA cuts Instance Quadratic model Epigraph model #nodes time gap #nodes time gap Fisher150iris 348,209 3600.00 1890.73 37080 3600.00 711.88 German22 7569 10.77 0.00 3247 12.81 0.00 German59 534,709 1949.13 0.00 90712 594.55 0.00 body-measurements 74,701 3600.00 15297.12 11605 3600.00 41450.30 cities-coord-202 218,931 3600.00 727.52 90,433 3600.00 187.54 cities-coord-666 26,602 3600.00 58814.95 21,894 3600.00 2596.30 concrete-compressive 1 3600.00 ∞1 3600.00 ∞ glass-identification 145,443 3600.00 1034.68 14,860 3600.00 947.23 image-segmentation 4374 3600.00 ∞304 3600.00 ∞ padberg-rinaldi-hole-dri 150,154 3600.00 ∞1027 3600.00 ∞ reinelt-hole-drilling 379,375 3600.00 ∞9789 3600.00 174671.81 ruspini 819,053 3600.00 871.99 216,584 2304.35 0.00 telugu-indian-vowel 281,326 3600.00 ∞9217 3600.00 97407.83 a1 111,630 3600.00 ∞558 3600.00 ∞ a2 30,873 3600.00 ∞747 3600.00 ∞ a3 0 3600.00 ∞533 3600.00 ∞ dim 7302 3600.00 ∞544 3600.00 ∞ g2-2-30 35,125 3600.00 99442.00 4829 3600.00 1247530.48 g2-2-50 34,063 3600.00 69060.41 5487 3600.00 2292990.33 g2-2-70 12,876 3600.00 ∞3349 3600.00 195477.72 s1 5624 3600.00 ∞1730 3600.00 ∞ s2 5496 3600.00 ∞837 3600.00 ∞ s3 1599 3600.00 ∞920 3600.00 ∞ s4 65,193 3600.00 ∞1004 3600.00 ∞ unbalance 5514 3600.00 ∞971 3600.00 ∞ 123 184 Journal of Global Optimization (2023) 87:133–189 Table 31 Results for 3 clusters using enabled heuristics, barycenter+convexity+cone+distance propagators, and OA cuts Instance Quadratic model Epigraph model #nodes time gap #nodes time gap Fisher150iris 446,923 3600.00 1849.96 59,165 3600.00 659.99 German22 7569 8.77 0.00 3247 10.52 0.00 German59 534,709 1597.28 0.00 91,402 503.05 0.00 body-measurements 83,144 3600.00 14592.81 15,451 3600.00 37113.59 cities-coord-202 239,312 3600.00 699.49 130,068 3600.00 147.72 cities-coord-666 27,921 3600.00 57591.83 29,912 3600.00 1974.12 concrete-compressive 1 3600.00 ∞1 3600.00 ∞ glass-identification 168,283 3600.00 1013.11 19,760 3600.00 817.90 image-segmentation 4830 3600.00 ∞393 3600.00 ∞ padberg-rinaldi-hole-dri 174,185 3600.00 ∞3265 3600.00 ∞ reinelt-hole-drilling 419,551 3600.00 ∞14,618 3600.00 52928.16 ruspini 906,635 3600.00 796.01 218,941 1775.63 0.00 telugu-indian-vowel 325,351 3600.00 ∞10,984 3600.00 97407.83 a1 143,309 3600.00 ∞699 3600.00 ∞ a2 41,088 3600.00 ∞803 3600.00 ∞ a3 0 3600.00 ∞1335 3600.00 ∞ dim 7947 3600.00 ∞682 3600.00 ∞ g2-2-30 35,530 3600.00 98431.46 5812 3600.00 1203432.46 g2-2-50 37,600 3600.00 63371.45 6492 3600.00 683571.57 g2-2-70 13,865 3600.00 ∞4075 3600.00 126844.00 s1 4586 3600.00 ∞1798 3600.00 ∞ s2 11,726 3600.00 ∞941 3600.00 ∞ s3 1599 3600.00 ∞997 3600.00 ∞ s4 52,430 3600.00 ∞1093 3600.00 ∞ unbalance 10,689 3600.00 ∞1340 3600.00 15265661.41 123 Journal of Global Optimization (2023) 87:133–189 185 Table 32 Results for 3 clusters using enabled heuristics, propagators, cuts, and entropy branching rule Instance quadratic model epigraph model #nodes time gap #nodes time gap Fisher150iris 343,703 3600.00 1477.66 71,751 3600.00 22683.70 German22 45,150 52.35 0.00 5333 12.16 0.00 German59 936,893 3600.00 56.04 336,433 1635.47 0.00 body-measurements 60,997 3600.00 14613.99 14,271 3600.00 16576.14 cities-coord-202 296,551 3600.00 301.34 155,242 3600.00 421.92 cities-coord-666 26,269 3600.00 7757.89 45,297 3600.00 2234.12 concrete-compressive 1 3600.00 ∞1 3600.00 ∞ glass-identification 169,038 3600.00 2547.47 19,694 3600.00 4728.44 image-segmentation 4722 3600.00 ∞1068 3600.00 ∞ padberg-rinaldi-hole-dri 171,210 3600.00 ∞11,198 3600.00 ∞ reinelt-hole-drilling 418,898 3600.00 ∞45,383 3600.00 ∞ ruspini 786,957 3600.00 384.48 393,872 3600.00 225.40 telugu-indian-vowel 317,636 3600.00 ∞8674 3600.00 194927.78 a1 138,799 3600.00 ∞1212 3600.00 ∞ a2 41,024 3600.00 ∞594 3600.00 ∞ a3 0 3600.00 ∞276 3600.00 ∞ dim 6776 3600.00 ∞1002 3600.00 ∞ g2-2-30 30,811 3600.00 72575.07 8337 3600.00 42157.54 g2-2-50 41,034 3600.00 61675.17 6234 3600.00 62328.59 g2-2-70 9287 3600.00 ∞5348 3600.00 174453.15 s1 4586 3600.00 ∞252 3600.00 ∞ s2 11,726 3600.00 ∞2999 3600.00 ∞ s3 1599 3600.00 ∞1978 3600.00 ∞ s4 52,430 3600.00 ∞3053 3600.00 ∞ unbalance 10,689 3600.00 ∞409 3600.00 ∞ 123 186 Journal of Global Optimization (2023) 87:133–189 Table 33 Results for 3 clusters using enabled heuristics, propagators, cuts, and distance branching rule Instance Quadratic model Epigraph model #nodes time gap #nodes time gap Fisher150iris 429,083 3600.00 2347.78 92,344 3600.00 1390.57 German22 11,515 15.42 0.00 2855 6.51 0.00 German59 295,515 1301.84 0.00 148,496 711.93 0.00 body-measurements 80,408 3600.00 29828.45 12,074 3600.00 ∞ cities-coord-202 254,608 3600.00 879.26 163,946 3600.00 470.39 cities-coord-666 39,469 3600.00 16406.07 74,088 3600.00 ∞ concrete-compressive 1 3600.00 ∞1 3600.00 ∞ glass-identification 167,142 3600.00 1986.31 17,442 3600.00 1406.33 image-segmentation 7073 3600.00 ∞3111 3600.00 ∞ padberg-rinaldi-hole-dri 175,960 3600.00 ∞10,875 3600.00 106261.12 reinelt-hole-drilling 411,784 3600.00 ∞5711 3600.00 ∞ ruspini 919,006 3600.00 578.31 221,022 1548.04 0.00 telugu-indian-vowel 340,173 3600.00 ∞15,303 3600.00 ∞ a1 141,729 3600.00 ∞7476 3600.00 ∞ a2 41,553 3600.00 ∞526 3600.00 ∞ a3 0 3600.00 ∞894 3600.00 ∞ dim 9647 3600.00 ∞1015 3600.00 ∞ g2-2-30 40,223 3600.00 231297.98 7569 3600.00 ∞ g2-2-50 37,147 3600.00 33690.28 5833 3600.00 ∞ g2-2-70 33,182 3600.00 ∞5532 3600.00 65019.56 s1 4586 3600.00 ∞4366 3600.00 ∞ s2 11,726 3600.00 ∞3905 3600.00 ∞ s3 1599 3600.00 ∞7583 3600.00 ∞ s4 52,430 3600.00 ∞10,129 3600.00 ∞ unbalance 10,398 3600.00 ∞1212 3600.00 ∞ References 1. Achterberg, T., Koch, T., Martin, A.: Branching rules revisited. Oper. Res. Lett. 33(1), 42–54 (2005). https://doi.org/10.1016/j.orl.2004.04.002 2. Aloise, D., Deshpande, A., Hansen, P., Popat, P.: NP-hardness of Euclidean sum-of-squares clustering. Mach. Learn. 75, 245–248 (2009). https://doi.org/10.1007/s10994-009-5103-0 3. Aloise, D., Hansen, P.: A branch-and-cut SDP-based algorithm for minimum sum-of-squares clustering. Pesquisa Operacional 29, 503–516 (2009). https://doi.org/10.1590/S0101-74382009000300002 4. Aloise, D., Hansen, P.: Evaluating a branch-and-bound RLT-based algorithm for minimum sum-of-squares clustering. J. Global Optim. 49, 449–465 (2011). https://doi.org/10.1007/s10898-010-9571-3 5. Aloise, D., Hansen, P., Liberti, L.: An improved column generation algorithm for minimum sum-ofsquares clustering. Math. Program. 131, 195–220 (2012). https://doi.org/10.1007/s10107-010-0349-7 6. Barber, C.B., Dobkin, D.P., Huhdanpaa, H.: The Quickhull algorithm for convex hulls. ACM Trans. Math. Softw. 22(4), 469–483 (1996). https://doi.org/10.1145/235815.235821 7. Brusco, M.J.: A Repetitive Branch-and-Bound Procedure for Minimum Within-Cluster Sums of Squares Partitioning. Psychometrika 71(2), 347–363 (2006). https://doi.org/10.1007/s11336-004-1218-1 8. Burgard, J.P., Costa, C.M., Schmidt, M.: Decomposition methods for Robustified k-means clustering problems: if less conservative does not mean less bad. Ann. Oper. Res. (2022). https://doi.org/10.1007/ s10479-022-04818-w 123 Journal of Global Optimization (2023) 87:133–189 187 9. Chen, C., Luo, J., Parker, K.: Image segmentation via adaptive Kmean clustering and knowledge-based morphological operations with biomedical applications. IEEE Trans. Image Process. 7(12), 1673–1683 (1998). https://doi.org/10.1109/83.730379 10. Cuesta-Albertos, J.A., Fraiman, R.: Impartial trimmed k-means for functional data. Comput. Stat. Data Anal. 51(10), 4864–4877 (2007). https://doi.org/10.1016/j.csda.2006.07.011 11. Dasgupta, S.: The hardness of k-means clustering. Tech. rep. Technical Report CS2008-0916. University of California, Department of Computer Science and Engineering. (2007). http://cseweb.ucsd.edu/~dasgupta/ papers/kmeans.pdf 12. Datta, S., Datta, S.: Comparisons and validation of statistical clustering techniques for microarray gene expression data. Bioinformatics 19(4), 459–466 (2003). https://doi.org/10.1093/bioinformatics/btg025 13. De Rosa, A., Khajavirad, A.: The ratio-cut polytope and K-means clustering. SIAM J. Optim. 32(1), 173–203 (2022). https://doi.org/10.1137/20M1348601 14. Deza, M.M., Laurent, M.: Geometry of Cuts and Metrics. Springer, Berlin (1997). https://doi.org/10. 1007/978-3-642-04295-9 15. Diehr, G.: Evaluation of a branch and bound algorithm for clustering. SIAM J. Sci. Stat. Comput. 6(2), 268–284 (1985). https://doi.org/10.1137/0906020 16. Dua, D., Graff, C.: UCI Machine Learning Repository. (2017). http://archive.ics.uci.edu/ml 17. Duran, M.A., Grossmann, I.E.: An outer-approximation algorithm for a class of mixed-integer nonlinear programs. Math. Program. 36(3), 307–339 (1986). https://doi.org/10.1007/BF02592064 18. du Merle, O., Hansen, P., Jaumard, B., Mladenovic, N.: An interior point algorithm for minimum sum-of-squares clustering. SIAM J. Sci. Comput. 21(4), 1485–1505 (1999). https://doi.org/10.1137/ S1064827597328327 19. Fisher, R.A.: The use of multiple measurements in taxonomic problems. Ann. Eugen. 7(2), 179–188 (1936). https://doi.org/10.1111/j.1469-1809.1936.tb02137.x 20. Fletcher, R., Leyffer, S.: Solving mixed integer nonlinear programs by outer approximation. Math. Program. 66(1), 327–349 (1994). https://doi.org/10.1007/BF01581153 21. Floudas, C., Aggarwal, A., Ciric, A.: Global optimum search for nonconvex NLP and MINLP problems. Comput. Chem. Eng. 13(10), 1117–1132 (1989). https://doi.org/10.1016/0098-1354(89)87016-4 22. Fränti, P., Sieranoja, S.: k-means properties on six clustering benchmark datasets. Appl. Intell. 48(12), 4743–4759 (2018). https://doi.org/10.1007/s10489-018-1238-7 23. Fränti, P., Sieranoja, S.: How much can k-means be improved by using better initialization and repeats? Pattern Recogn. 93, 95–112 (2019). https://doi.org/10.1016/j.patcog.2019.04.014 24. Fukuda, K.: cdd/cdd+ Reference Manual. In: Institute for Operations Research, ETH-Zentrum, pp. 91–111 (1997) 25. Fukunaga, K., Narendra, P., Koontz, W.: A branch and bound clustering algorithm. IEEE Trans. Comput. 24(09), 908–915 (1975). https://doi.org/10.1109/T-C.1975.224336 26. Gamrath, G., Anderson, D., Bestuzheva, K., Chen, W.-K., Eifler, L., Gasse, M., Gemander, P., Gleixner, A., Gottwald, L., Halbig, K., Hendel, G., Hojny, C., Koch, T., Le Bodic, P., Maher, S.J., Matter, F., Miltenberger, M., Mühmer, E., Müller, B., Pfetsch, M.E., Schlösser, F., Serrano, F., Shinano, Y., Tawfik, C., Vigerske, S., Wegscheider, F., Weninger, D., Witzig, J.: The SCIP Optimization Suite 7.0. eng. Tech. rep. 20-10. Takustr. 7, 14195 Berlin: ZIB (2020) 27. Gilpin, A., Sandholm, T.: Information-theoretic approaches to branching in search. Discrete Optim. 8(2), 147–159 (2011). https://doi.org/10.1016/j.disopt.2010.07.001 28. Gonzalez, T.F.: Clustering to minimize the maximum intercluster distance. Theor. Comput. Sci. 38, 293– 306 (1985). https://doi.org/10.1016/0304-3975(85)90224-5 29. Grötschel, M.H.: Solution of large-scale symmetric travelling salesman problems. Math. Program. 51, 141–202 (1991). https://doi.org/10.1007/BF01586932 30. Guns, T., Dao, T.-B.-H., Vrain, C., Duong, K.-C.: Repetitive branch-andbound using constraint programming for constrained minimum sum-of-squares clustering. In: Proceedings of the Twenty-second European Conference on Artificial Intelligence (ECAI’16). IOS Press, NLD, pp. 462–470 (2016). https:// doi.org/10.3233/978-1-61499-672-9-462 31. Han, S.: Spatial stratification and socio-spatial inequalities: the case of Seoul and Busan in South Korea. Human. Soc. Sci. Commun. 9(1), 23 (2022). https://doi.org/10.1057/s41599-022-01035-5 32. He, H., Chen, J., Jin, H., Chen, S.-H.: Trading strategies based on K-means clustering and regression models. In: Chen, S.-H., Wang, P.P., Kuo, T.-W. (eds.), Computational Intelligence in Economics and Finance: Volume II, pp. 123–134. Springer, Berlin (2007). https://doi.org/10.1007/978-3-540-728214_7 33. Heinz, G., Peterson, L.J., Johnson, R.W., Kerk, C.J.: Exploring relationships in body dimensions. J. Stat. Educ. (2003). https://doi.org/10.1080/10691898.2003.11910711 123 188 Journal of Global Optimization (2023) 87:133–189 34. Horst, R., Tuy, H.: Global Optimization. Springer, Berlin (1996). https://doi.org/10.1007/978-3-66203199-5 35. Hua, K., Shi, M., Cao, Y.: A Scalable deterministic global optimization algorithm for clustering problems. In: International Conference on Machine Learning. PMLR, pp. 4391–4401 (2021). https://proceedings. mlr.press/v139/hua21a.html 36. Kaibel, V., Peinhardt, M., Pfetsch, M.E.: Orbitopal fixing. Discret. Optim. 8(4), 595–610 (2011). https:// doi.org/10.1016/j.disopt.2011.07.001 37. Kaibel, V., Pfetsch, M.E.: Packing and partitioning orbitopes. Math. Program. 114(1), 1–36 (2008). https:// doi.org/10.1007/s10107-006-0081-5 38. Liberti, L., Manca, B.: Side-constrained minimum sum-of-squares clustering: mathematical programming and random projections. J. Global Optim. (2021). https://doi.org/10.1007/s10898-021-01047-6 39. Lloyd, S.: Least squares quantization in PCM. IEEE Trans. Inf. Theory 28(2), 129–137 (1982). https:// doi.org/10.1109/TIT.1982.1056489 40. MacQueen, J.: Some methods for classification and analysis of multivariate observations. In: Proceedings of the Fifth Berkeley Symposium on Mathematical Statistics and Probability, Volume 1: Statistics, pp. 281– 297. University of California Press, Berkeley (1967). https://projecteuclid.org/euclid.bsmsp/1200512992 41. Mahajan, M., Nimbhorkar, P., Varadarajan, K.: The planar k-means problem is NP-hard. In: Theoretical Computer Science 442. Special Issue on the Workshop on Algorithms and Computation (WALCOM 2009), pp. 13–21 (2012). https://doi.org/10.1016/j.tcs.2010.05.034 42. Padberg, M., Rinaldi, G.: A branch-and-cut algorithm for the resolution of large-scale symmetric traveling salesman problems. SIAM Rev. 33(1), 60–100 (1991). https://doi.org/10.1137/1033004 43. Pal, S.K., Majumder, D.D.: Fuzzy sets and decision making approaches in vowel and speaker recognition. IEEE Trans. Syst. Man Cybern. 7(8), 625–629 (1977). https://doi.org/10.1109/TSMC.1977.4309789 44. Peng, J., Wei, Y.: Approximating k-means-type clustering via semidefinite programming. SIAM J. Optim. 18(1), 186–205 (2007). https://doi.org/10.1137/050641983 45. Peng, J., Xia, Y.: A cutting algorithm for the minimum sum-of-squared error clustering. In: Proceedings of the 2005 SIAM International Conference on Data Mining, pp. 150–160 (2005). https://doi.org/10.1137/ 1.9781611972757.14 46. Peng, J., Xia, Y.: A new theoretical framework for k-means-type clustering. In: Foundations and advances in data mining. Springer, Berlin, pp. 79–96 (2005). https://doi.org/10.1007/11362197_4 47. Piccialli, V., Sudoso, A.M., Wiegele, A.: SOS-SDP: an exact solver for minimum sum-of-squares clustering. INFORMS J. Comput. 34(4), 2144–2162 (2022). https://doi.org/10.1287/ijoc.2022.1166 48. Plastria, F.: Formulating logical implications in combinatorial optimisation. Eur. J. Oper. Res. 140(2), 338–353 (2002). https://doi.org/10.1016/S0377-2217(02)00073-5 49. Prasad, M.N., Hanasusanto, G.A.: Improved conic reformulations for k-means clustering. SIAM J. Optim. 28(4), 3105–3126 (2018). https://doi.org/10.1137/17M1135724 50. Quesada, I., Grossmann, I.E.: An LP/NLP based branch and bound algorithm for convex MINLP optimization problems. Comput. Chem. Eng. 16(10–11), 937–947 (1992). https://doi.org/10.1016/00981354(92)80028-8 51. Reinelt, G.: TSPLIB-A traveling salesman problem library. ORSA J. Comput. 3(4), 376–384 (1991). https://doi.org/10.1287/ijoc.3.4.376 52. Ruspini, E.H.: Numerical methods for fuzzy clustering. Inf. Sci. 2(3), 319–350 (1970). https://doi.org/ 10.1016/S0020-0255(70)80056-1 53. Sangalli, L.M., Secchi, P., Vantini, S., Vitelli, V.: k-mean alignment for curve clustering. Comput. Stat. Data Anal. 54(5), 1219–1233 (2010). https://doi.org/10.1016/j.csda.2009.12.008 54. Shannon, C.E.: A mathematical theory of communication. Bell Syst. Tech. J. 27(3), 379–423 (1948). https://doi.org/10.1002/j.1538-7305.1948.tb01338.x 55. Sherali, H.D., Desai, J.: A global optimization RLT-based approach for solving the hard clustering problem. J. Global Optim. 32, 281–306 (2005). https://doi.org/10.1007/s10898-004-2706-7 56. Sobol’, I.: On the distribution of points in a cube and the approximate evaluation of integrals. USSR Comput. Math. Math. Phys. 7(4), 86–112 (1967). https://doi.org/10.1016/0041-5553(67)90144-9 57. Späth, H.: Cluster Analysis Algorithms for Data Reduction and Classification of Objects. Horwood, Bristol (1980) 58. Steinley, D.: K-means clustering: a half-century synthesis. Br. J. Math. Stat. Psychol. 59(1), 1–34 (2006). https://doi.org/10.1348/000711005X48266 59. Tan, M.P., Broach, J.R., Floudas, C.A.: A novel clustering approach and prediction of optimal number of clusters: global optimum search with enhanced positioning. J. Global Optim. 39, 323–346 (2007). https:// doi.org/10.1007/s10898-007-9140-6 60. Tïrn˘auc˘a, C., Gómez-Pérez, D., Balcázar, J.L., Montaña, J.L.: Global optimality in k-means clustering. Inf. Sci. 439–440, 79–94 (2018). https://doi.org/10.1016/j.ins.2018.02.001 123