Mixed-Integer Programming Techniques for the Minimum Sum-of-Squares Clustering Problem

. 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 diﬀerent 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 signiﬁcantly 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.


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 (Cuesta-Albertos and Fraiman 2007; Sangalli et al. 2010), image processing (Chen et al. 1998), bio-informatics (Datta and Datta 2003), economics (He et al. 2007), and social sciences (Han 2022).For a detailed survey of the history of clustering problems we refer to Steinley (2006).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. (2009), Dasgupta (2007), and Mahajan et al. (2012).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 (1982) and MacQueen (1967).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-andbound methods is presented by Fukunaga et al. (1975), which has been refined later on by Diehr (1985).A variant of a so-called repetitive branch-and-bound method has been devised by Brusco (2006), 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 (2005).The authors use reformulationlinearization-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 (2002), which is, generally speaking, a modeling tutorial paper but which also contains a discussion of symmetry breaking constraints for the clustering problem.Aloise and Hansen (2011) also consider the MSSC problem and try to re-produce the results of Sherali and Desai (2005).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. (1999) and which has been re-considered and improved by Aloise et al. (2012).Also other classic techniques of mixed-integer (non)linear optimization have been applied such as generalized Benders decomposition in Floudas et al. (1989) and Tan et al. (2007).Alternatively, Peng and Xia (2005a) consider the MSSC problem as a concave minimization problem and adapt Tuy's cut method, see Horst and Tuy (1996), for solving the problem.Further, Prasad and Hanasusanto (2018) propose improved conic reformulations of the MSSC problem and also study some symmetry breaking techniques.Tîrnăucă et al. (2018) 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 (2007) and Peng and Xia (2005b) proved the equivalence between the MSSC problem and a 0-1 SDP reformulation.Based on this 0-1 SDP model, Aloise and Hansen (2009) propose a branch-and-cut algorithm and solve instances with up to 202 data points to global optimality.More recently, Piccialli et al. (2022) 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-ofthe-art branch-and-bound algorithm for the MSSC problem.SDP-like models have also been used in De Rosa and Khajavirad (2022) to use Z = XX ∈ [0, 1] n×n for 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 (1997) 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. (2021) and Liberti and Manca (2021).Besides that, in Liberti and Manca (2021), 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 (Section 2), develop cutting planes (Section 3), propagation methods (Section 4), as well as problem-specific branching rules (Section 5) and primal heuristics (Section 6).We implement and test all techniques in the open-source MINLP solver SCIP; see Gamrath et al. (2020).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 Section 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 Section 9. Our code is publicly available at GitHub1 .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-of-the-art and SDP-based techniques as studied in Piccialli et al. (2022).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.

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 ⊆ R d and 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 ∈ P to 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 c j is minimal.This problem is modeled via the following MINLP: p∈P j∈ [k] x pj p − c j 2 (1a) j∈ [k] x pj = 1, p ∈ P, (1b) The binary variables x pj are the assignment variables that model whether the data point p is assigned to cluster j (x pj = 1) or not (x pj = 0).Moreover, B ⊆ R d is a set that contains all points in P .This can, e.g., be the bounding box of } 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 x with 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 j∈ [k] x pj = 1, p ∈ P, (2c) where M p are sufficiently large numbers.The objective function is linear now and we obtain the additional quadratic and convex constraints in (2b).
For every p ∈ P , M p can be chosen to be the maximum distance of p to any other point p ∈ P .An overestimation can be easily computed via where i and u i are 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 c j , j ∈ [k], is p∈P p x pj p∈P x pj , i.e., the barycenter of all points assigned to cluster j.

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. (2012) as well as Sherali and Desai (2005).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): 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 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 κ j with range {k − 1, . . ., |P | − 1} and link it with the x-variables via the linear constraint That is, κ j describes the number of data points that are not assigned to cluster j.If a lower bound κ j and an upper bound κj on κ j is given, this equation implies the inequalities That is, they describe localized versions of cardinality cuts that get stronger if xvariables get fixed.Another side effect of the auxiliary variables κ j is 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): This allows to solve Model (2) in an outer approximation or LP/NLP-based branchand-bound fashion; see Duran and Grossmann (1986) and Fletcher and Leyffer (1994) or Quesada and Grossmann (1992), 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 (η, ĉ, x) of Model (2).In the original outer approximation method, this is done by fixing the binary variables x in Model (2) and solving the resulting convex NLP subproblem; see Duran and Grossmann (1986).
The benefit in our specific application is that solving the subproblem boils down to a simple computation of the barycenters ĉ, 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 (η, ĉ, x) of the subproblem, it holds ηpj for all feasible points (η, c, x) of the updated relaxation.In other words, adding the outer-approximation cuts bounds the optimal objective value of the relaxation with fixed binaries x = x from below by p∈P j∈[k] ηpj .Consequently, the updated relaxation yields a solution with a new, previously unseen, cluster assignment x or the optimality gap is closed.Thus, iterating this process terminates after a finite number of steps; see Duran and Grossmann (1986) or Fletcher and Leyffer (1994) 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.

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], we denote by P j ⊆ P the set of all data points p ∈ P whose corresponding variable x pj has been fixed to 1 at the current node of the branchand-bound tree.That is, we have already decided to assign p to cluster j.Moreover, we denote by P j ⊆ P all data points p such that x pj has not been fixed to 0 yet, i.e., p is already or can still be assigned to cluster j.Note that P j ⊆ P j .For a continuous variable z, i.e., for the c-and η-variables, we denote by z and z the lower and upper bound on z at the current node, respectively.4.1.Barycenter Propagation.Given a non-empty set of data points Q ⊆ P defining a cluster, the optimal choice for its centroid is the barycenter 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 Consequently, a lower bound on the objective is given by j∈[k] D(P j ).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 (2006), see also Guns et al. (2016), 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(P j ).Otherwise, if p / ∈ P j and η 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.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 c j i with (i, j) , we compute γ j,r = C(P j ∪ {p 1 , . . ., p r }), i.e., the barycenter of the data points contained in P j and the data points with the r smallest ith coordinates that are not contained in P j .As we show next, the ith coordinates of these barycenters can be used to compute a lower bound on c j i .Lemma 4.1.A valid lower bound on c j i is given by min r∈[s]0 γ j,r i .Proof.Let Q ⊆ P j \P j and assume |Q| = r.Then, C(P j ∪Q) i ≥ C(P j ∪{p 1 , . . ., p r }) i , because the points p 1 , . . ., p r are points with the r smallest ith coordinates.Consequently, to find a lower bound on the centroids, it is sufficient to consider C(P j ∪ {p 1 , . . ., p r }) for each r ∈ [s] 0 .
Analogously, an upper bound is given by max r∈[s]0 C(P j ∪ {p s−r , . . ., p s }) i .Since computing an iterative sequence of barycenters can be done using the formula we can compute the minimum and maximum values for all coordinates and clusters in O(kd|P |) time.

Convexity and Cone Propagation.
Based on optimality arguments, we can also derive rules to assign data points p ∈ P j \ P j to 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 P 1 , . . ., P k such that, for each j ∈ [k], we have conv(P j ) ∩ P = P j .
Proof.Given an optimal allocation of the k centroids, the Voronoi cells cover the entire R d and 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 P 1 ⊆ C 1 in any optimal solution.If there exist p ∈ P \ P 1 that are contained in C 1 , they are necessarily contained in the boundary of C 1 .Hence, if we change the assignment of these points to P 1 , this does not change the objective of MSSC.The assertion thus holds for P 1 , and we can use the same arguments iteratively to conclude the proof.
As a consequence, the convexity propagator computes conv(P j ) for each j ∈ [k].If there exists p ∈ P ∩ conv(P j ) it performs the following steps: If p / ∈ P j holds, 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, x pj 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 Figure 1.q Figure 1.Illustration of cone-based propagation.If q is not contained in the red cluster, none of the black points can be contained in the red cluster.Assigning the white points to the red cluster is still possible.
Lemma 4.3.Let P 1 ∪ • • • ∪ P k be a partition of a finite set P ⊆ R d .Suppose conv(P j ) ∩ P = P j for each j ∈ [k].Then, for every q ∈ P \ P j , q + cone{−(p − q) : p ∈ P j } ⊇ P \ P j .
Proof.Note that q + cone{p − q : p ∈ P j } is the smallest cone with apex q that contains conv(P j ), because q / ∈ P j and we shoot rays from q through each of the finitely many points in P j .Hence, negating these rays leads to a cone that cannot contain any point from P j as it does not contain conv(P j ).
We can use this observation as follows.If there is q ∈ P \ P j , i.e., x qj is fixed to 0, then all points in q + cone{−(p − q) : p ∈ P j } can be fixed to 0 as well.
In arbitrary dimensions, the convexity propagator cannot be implemented efficiently, because conv(P j ) might have Ω(2 d ) many facets.In small dimensions, computing convex hulls can be done rather quickly and, as our numerical results will indicate, have a very positive impact on the time needed to solve MSSC.

Distance Propagation.
The distance propagator provides another set of rules to fix variables x pj , (p, j) ∈ P × [k], to 0. To this end, it defines for each j ∈ [k] the bounding box } for the centroids.That is, the smallest box that contains the centroid for cluster j based on local variable bound information.Afterward, for each p ∈ P and j ∈ [k], it computes the minimum and maximum distances D min j,p and D max j,p to the bounding box B j , i.e., D min j,p = min { p − x : x ∈ B j } , D max j,p = max { p − x : x ∈ B j } .Since a data point is assigned to a centroid of minimum distance in an optimal solution, p cannot be assigned to cluster j ∈ [k] if there is j ∈ [k] with D max j ,p < D min j,p .Consequently, x pj can be fixed to 0 in this case.
Finding B j and computing the minimum and maximum distance of a point to a box can be done in O(d) time.Hence, the distance propagator runs in O(kd|P |) time.

Branching rules
After a node of the branch-and-bound tree has been processed (adding cutting planes, propagation), branching rules split the current subproblem into further subproblems to, e.g., tighten the problem formulation or enforce integrality of variables.The decision on how to split the current subproblem is typically guided by the solution of the subproblem's LP relaxation.To enforce integrality of variables (in the notation of the MSSC problem), one typically selects a variable x pj whose value in the LP solution is non-integral.Then, two subproblems are created that fix x pj to 0 and 1, respectively.Despite the existence of many branching rules that perform well for generic problems, see Achterberg et al. (2005), there also exist branching rules tailored to a specific problem.This becomes relevant, because they might allow to derive further reductions based on problem structure or further components of a solver such as propagation mechanisms.
For integer programs, Gilpin and Sandholm (2011) proposed four families of branching rules motivated from an information-theoretic perspective.The common ground of all their rules is to interpret the values x pj as the probability that a point p ∈ P is assigned to cluster j ∈ [k].Using their rules, they aim at reducing the assignment uncertainty in the current subtree.The first and second rule use a look-ahead approach, similar to strong branching.For the MSSC problem with a large number of clusters or points, look-ahead branching rules may become computationally prohibitive when applied to the full problem.Hence, we do not use them.The third family is called entropic look-ahead-free variable selection and the fourth is its extension for a multi-variable branching version.Since we think that the third family might be helpful for the MSSC problem, we describe it in the following.Afterward, two novel branching rules for the MSSC problem will be presented.
5.1.Entropy Branching.Suppose that the optimal LP solution of the current subproblem does not satisfy the integrality constraints, which means that the relaxed solution x is non-integral.Let X be the set of all branching candidates that are non-integral in this LP solution, i.e., Since each of these xpj is non-integral, the cluster assignment of point p is not fixed.
For the assignment to be fixed, xpj has to be one.Due to Constraints (1b) or (2c), xpj can be seen as a kind of posterior probability of point p to belong to cluster j (Zheng et al. 2017).A good strategy for branching would be to select the point p for which the probabilities for each cluster assignment are almost the same.
The most unclear situation is where xpj = 1 k holds for all j ∈ [k].Here, each cluster assignment is equally probable for point p.This can be seen as a homogeneous information setting.The level of homogeneity can be measured via the Shannon entropy of point p (Shannon 1948).More precisely, for each variable xpj ∈ X the entropy of point p with probabilities xpj , j ∈ [k], is The maximum entropy (H p = log 2 k) occurs in the above mentioned extreme case.That is, the current LP solution does not provide any information on the best (or most probable) cluster assignment of point p.The minimum entropy is obtained when there is a clear cluster assignment.In this situation, let the point p be already assigned to a cluster, e.g., to cluster j = 1 and hence xp1 = 1.Due to (1b) (or (2c)), xpj = 0 for j ∈ [k] \ {1}.The entropy of p is then where 0 log 2 0 is taken to be zero.
We are interested in finding the point corresponding to the fractional variable xp * j * ∈ X such that the entropy of point p * is the maximal over all points with fractional variables, i.e., we search for the most uncertain assignment.The point is formally given by p * ∈ arg max For this point p * , we select the cluster index j to branch on arbitrarily, i.e., we create two subproblems by adding either xp * j = 1 or xp * j = 0. 5.2.Distance Branching.Next, we describe three branching rules with a geometric motivation.The first one, called distance branching, is based on the intuition that clusters should be rather compact (opposed to being spread out).Given the current LP solution with its suggestion for the centroids c j , j ∈ [k], the variable x pj ∈ X selected for branching is the one corresponding to the data point p and cluster j that are most apart from each other, i.e., we find Then, we branch on the fractional variable xp * j * , creating two subproblems by adding either xp * j * = 1 or xp * j * = 0.If an optimal cluster is indeed compact, then the 0-subproblem contains an optimal solution.Otherwise, the convexity propagator has the potential to also fix additional variables that lie between p * and the remaining points of cluster j * in the 1-subproblem.5.3.Centrality Branching.Since the distance branching rule is tailored towards the extremes of compact vs. far spread-out clusters, the centrality branching rule takes a more balanced approach by selecting a point whose distance to a cluster is not too big.Given the current LP solution with its suggestion for the centroids c j , j ∈ [k], we would like to branch on the non-integral variable x pj corresponding to the data point p and cluster j that is lying in the center of the cloud of unassigned data points.To obtain a cheap evaluation, we take the point p * that is in the center of all centroids, i.e., From this point p * , we select an arbitrary variable xp * j * , creating two subproblems by adding either xp * j * = 1 or xp * j * = 0.If the distance of p * to cluster j * is not too small, then there is a chance that the convexity propagator can fix further data points to be contained in cluster j * in the 1-subproblem.Opposed to the distance branching rule, however, also the 0-subproblem becomes relevant for the convexity propagator as it might fix some variables to 0 based on its cone propagation.
5.4.Pairs-in-the-Middle-of-Pairs Branching.The next branching rule is a variation of the last one, but now we branch on more general linear inequalities.
Given the current LP solution with its suggestion for the centroids c j , j ∈ [k], we would like to branch on the sum of a pair of non-integral variables corresponding to the two data points located in the middle of a pair of clusters.First assume that there are only two clusters with corresponding centroids c 1 and c 2 .We want to find the points that are nearest to the point lying half way from centroid c 1 to centroid c 2 .Any point p ∈ P that lies on the line segment between c 1 and c 2 , minimizes the sum p − c 1 + p − c 2 by the triangle inequality of the Euclidean distance.With the same reasoning, the smaller this sum is, the nearer is point p to the line segment between the two clusters.However, there can be multiple points p with the same value for p − c 1 + p − c 2 .We are interested in the ones that are nearest to the middle.Thus, we penalize longer distances by minimizing the sum p − c 1 2 + p − c 2 2 of squares instead.Now assume that there are more than two clusters which have pairwise the exact same distance of centroids to each other.Then, still looking for the two points p ∈ P that minimize p − c j 2 + p − c j 2 for j, j ∈ [k] with j = j gives us the desired points.Let p and q be the selected points and j and j be the selected clusters.We then compute the sums in the current LP solution, xpj + xqj and xpj + xqj , and select the sum that is most fractional or least fractional.Both versions are tested in our numerical experiments.Suppose that the first sum is selected, which means that the selected cluster is j.Then, we branch on

Primal Heuristics
Primal heuristics try to find feasible solutions of good quality in a short amount of time.Having good feasible solutions at hand early in the solving process is crucial.Feasible solutions help to prune branch-and-bound nodes based on bounding as well as to perform further fixings and reductions.Moreover, a user may already be satisfied with the quality of the heuristic solution, such that the solving process can be stopped at an early stage.In this section we present three primal heuristics for the MSSC problem.
6.1.A Root-Node Heuristic.To obtain a first feasible point, i.e., a point for warm-starting, we use the k-means algorithm, which is the most popular heuristic for finding a feasible solution for the MSSC problem; see, e.g., Lloyd (1982) and MacQueen (1967).It consists of two main steps.First, given an initial guess for the location of the centroids, each data point is assigned to the nearest centroid.Afterward, each centroid is updated by calculating the mean of the data points assigned to this centroid.This process is repeated until the centroids no longer change.To obtain an initial guess for the location of the centroids, we use the "furthest point heuristic", also known as "Maxmin" (Gonzalez 1985).The idea is to select the first centroid randomly within the respective bounding box and then obtain new centroids one by one.In each iteration, the next centroid is the point that is the furthest (max) from its nearest (min) existing centroid.Here, we choose the first data point as the first centroid.For a comparison of several initialization heuristics, see, e.g., Fränti and Sieranoja (2019).

6.2.
A Rounding Heuristic.Feasible solutions can be obtained at each node by applying a rounding scheme to the LP solution.We use the rounding heuristic proposed by Sherali and Desai (2005).For completeness, we also describe it here.
Given a non-integral LP solution (x, c), or (x, c, η) for Model (2), at a node of the branch-and-bound tree, we round the non-integral x-solution to the closest feasible binary solution x, while respecting the decisions that have already been made, i.e., if a data point is already assigned to a cluster it will remain in that cluster.First, we ensure that there are no empty clusters by finding, for each j ∈ [k] with P j = ∅, the point p ∈ P \ ∪ j∈[k] P j such that p ∈ arg max xpj : p ∈ P \ ∪ j∈[k] P j , and setting xpj = 1.To break a tie, the point with smallest index is chosen.Now, to ensure that the point p is only in one cluster, we set xpj = 0, for all j ∈ [k] \ {j}.
Furthermore, for each data point is not yet rounded, we find a cluster j * such that xpj * ∈ max {x pj : j ∈ [k]}.Again, we break ties by selecting the cluster with the smallest index.Then, we set xpj * = 1 and xpj = 0 for all j ∈ [k] \ {j * }.With x at hand, we can then compute the centroids for each cluster, see Observation 2.1, and obtain a feasible solution for the MSSC problem.6.3.An Improvement Heuristic.Given a feasible solution (x, c) of Model (1) or (x, c, η) of Model (2), we try to improve this solution by evaluating the loss function (i.e., the intra-variance) within each cluster.For that, consider the weighted value of the loss function restricted to cluster C j as It may happen that some clusters have a large loss function value, while some other clusters may have a very small one.Thus, we may find a better solution-regarding the sum of all losses-by splitting a cluster into two smaller clusters and joining two other clusters.This heuristic has been proposed in Burgard et al. (2022), where the motivation for its development is explained in more details.
The procedure is described as follows.For each pair of clusters (C j1 , C j2 ), we compute their joint centroid and the corresponding total loss via

Now, consider the set
which is the set of all possible combinations of three clusters such that the total loss within two joined clusters is smaller than the total loss within a third cluster.Note that the set Ψ can be empty.If so, this means that we cannot obtain a better solution by joining two clusters and splitting another one.On the other hand, i.e., if there exists (C j1 , C j2 , C j3 ) ∈ Ψ, then the total loss of the joined clusters C j1 and C j2 is smaller than the total loss within cluster C j3 .Thus, we obtain a better solution by joining C j1 and C j2 and by splitting cluster C j3 into two smaller clusters.To this end, we update the centroids in such a way that the clusters C j1 and C j2 are now one cluster with centroid c j1j2 , cluster C j3 receives two new centroids, and the other centroids remain the same, i.e., where c and c are obtained as follows.First we find the two furthest points in C j3 to be the initial guesses for the location of the centroids, i.e., Next, each point in C j3 is assigned to the closest centroid, either c or c .Then, the centroids c and c are updated based on this assignment.Now, the update of the assignments and centroids is repeated until they do not change anymore.This way we obtain the new centroids c and c that give us the desired splitting of cluster C j3 .Finally, if the set Ψ has more than one element, then we repeat the process starting with the element (C j1 , C j2 , C j3 ) that gives the minimum ratio F j1j2 /F j3 .Each time an element (C j1 , C j2 , C j3 ) is used, we exclude all the elements that contain C j1 , C j2 , or C j3 , because these clusters have already been modified.
With ĉ at hand, we can easily compute x and, thus, a new feasible solution (x, ĉ) is obtained.If the objective function value is better at this new solution, then we have found an improved solution out of (x, c).

Symmetry breaking
Note that both Model (1) and ( 2) are symmetric with respect to cluster assignments.That is, once a feasible solution has been found, one can generate equivalent (symmetric) solutions by exchanging the labels of the clusters.Such symmetries are known to deteriorate the performance of search-based approaches like branch-andbound, because symmetric subproblems are created repeatedly without providing the solver with new information.Such cluster symmetries can be handled in both models by imposing additional restrictions on the x-variables.If we interpret x as a binary matrix whose columns are labeled by clusters, then we can handle symmetries by enforcing that the columns of x need to be sorted lexicographically non-increasing.Since each row of matrix x has exactly one 1-entry due to (1b) and ( 2c), the lexicographic sorting can be imposed by orbitopal fixing, see Kaibel et al. (2011), and separating the symmetry handling inequalities developed by Kaibel and Pfetsch (2008).

Numerical Experiments
In this section, we report extensive computational results that show the benefits of the techniques proposed in Sections 3-6.To this end, we have incorporated all these techniques into the state-of-the-art solver SCIP; see Gamrath et al. (2020).As a reference for comparison, we use plain SCIP for both problem formulations (1) and (2).That is, we solve the two formulations without our problem-specific enhancements but with enabled symmetry handling.
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 setup are described in Sections 8.1 and 8.2, respectively.Then, in Section 8.3, we start the discussion 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 Section 8.4, we repeat the discussion but for the case k = 3.Finally, in Section 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 (2011) and Sherali and Desai (2005), as well as in Aloise et al. (2012).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 (2018).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 [−10 3 , 10 3 ] d .More precisely, for each dimension i ∈ [d], we compute the maximum and minimum coordinate value obtaining vi and v i , respectively.Then, we take If we do this for all data points, they get centered around the origin.Now, if w i = v i − u i < −10 3 or wi = vi − u i > 3 holds, we re-scale the data.The desired new bounds then are z i = −10 3 and zi = 10 3 .Thus, the re-scaled data point p is The corresponding instances that needed to be re-scaled are s1, s2, s3, s4, and unbalance.
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 Sections 3-6 are implemented as SCIP plugins written in C/C++ and our code is publicly available at GitHub 2 (githash 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 Section 7. To compute convex hulls and cones in the convexity propagator proposed in Section 4.2, we use the Qhull 3 C++ interface proposed by Barber et al. (1996).
We have also conducted experiments using the CDD library (Fukuda 1997) 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 1 h per instance.
2 https://github.com/christopherhojny/globally-solving-MSSC 3http://www.qhull.org.In the following, we discuss the impact of our techniques on solving the MSSC problem for k = 2 and k = 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 3600 s 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 t 1 , . . ., t n are provided as shifted geometric means n i=1 (t i + s) 1 /n − s to 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.
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.
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 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 3 corresponds 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.
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 Figure 2, where we plot the instances vs. the corresponding gap between the primal and dual bounds.This already demonstrates the great benefit that the  .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 10 4 , then it assumes the gap limit of 10 4 in the plot.
barycenter propagator adds to the solution process.As discussed in Section 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 Figure 2 again, we also see that for many instances the gaps improve.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 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.
Barycenter+Convexity+Cone Propagators.Although the barycenter and convexitycone 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 Figure 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.
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 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 Figure 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.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.
Branching Rules.The last components to be tested are branching rules.We have implemented all branching rules described in Section 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-themiddle-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" we refer to SCIP'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.
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 Figure 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 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.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.
Primal Heuristics.The summary results of plain SCIP and SCIP enabled with our primal heuristics are presented in Table 7.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 Figure 6.For this reason, we enable heuristics in the following experiments.
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 = 3 is much harder to solve than for k = 2, the propagation techniques are still helpful; see also Figure 7 where we compare the gaps obtained by using propagators.Therefore, in what follows, we discuss the benefits of the separate propagators in turn.Barycenter Propagator.By enabling only the barycenter propagator, we see that more nodes can be explored in the branch-and-bound tree.Particularly, when the quadratic model is used, there is a huge difference in the number of nodes when comparing the first two rows of Table 8.Moreover, the barycenter propagator allows one more instance to be solved, which could not be solved before by the quadratic model, resulting in a much lower mean running time.Finally, although the number of solved instances increases only slightly, we can see that the barycenter propagator has also a very positive effect on the other instances: the primal-dual gap improves significantly, in particular, for the quadratic model; see Figure 7.
Besides the improvement in the gaps for the epigraph model, the barycenter propagator also allows more nodes to be explored while reducing the running times considerably.Therefore, also for k = 3, the barycenter propagator allows to simplify the relaxations used in branch-and-bound.
Convexity+Cone Propagator.If we use only the convexity propagator and compare the results with plain SCIP and enabled heuristics, then we see that more nodes can be processed in significantly less time.Additionally enabling the cone propagator also allows to solve two more instances in the epigraph model; see Figure 7 again.Regarding the quadratic model, using the convexity and cone propagators allows to process more nodes in the branch-and-bound tree.This, however, comes at the price that the solvable instance requires more time.We conclude that the convexity and cone propagators enhance the solution process, where the effect is more dominant for the epigraph model.
Barycenter+Convexity+Cone Propagator.The quadratic model clearly benefits from using the barycenter propagator together with convexity and cone propagation as one more instance can be solved and in less time.For the epigraph model, no significant change regarding running times can be observed, however it has a positive effect on some gaps; see Figure 7.
Cutting Planes.Regarding the cutting planes, we again focus only on the OA cuts, since the localized version of the cardinality cuts does not positively affect the solution process in general.In Table 9 we present the aggregated results for the epigraph model with and without OA cuts.Note that by using the OA cuts we can solve the same number of instances, while requiring less time and much less nodes.The OA cuts help mainly in terms of dual bounds; see also Figure 8.Therefore, we conclude that the OA cuts are also beneficial for the epigraph model in the harder case of k = 3.   Distance Propagator.The distance propagator does not impact the solution process if used standalone.The reason is that, for most of the instances, the propagator does not find any reduction, which is to be expected because fixing x-variables to 0 is less powerful than fixing them to 1 (as the convexity propagator does, for instance).However, as more reductions and fixings are being made by other components, the distance propagator slightly helps; see the summary results shown in Table 10.
The distance propagator improves the solution process of both models, since more nodes can be explored and in less time while solving the same number of instances.Therefore, we enable this propagator as well in the following experiments.
Branching Rules.Finally, we report on the effect of branching rules.Again, the centrality and pairs-in-the-middle-of-pairs rule hinder the solution process.Therefore, we focus only on the entropy and distance branching rules, which are summarized in Table 11.The branching rules do not really change the results.In Figure 9, we see that the entropy branching rule has a positive effect on four gaps of the quadratic model, but harms the solution process of two other instances.Moreover, it requires more time.For the epigraph model, it is visibly bad.The standard and distance branching rules perform equally good for the quadratic model; see Figure 9 again.While, for the epigraph model, the standard branching rule performs best.Best Setting.To finalize the discussion of the numerical results for the case k = 3, we again compare plain SCIP with the best setting so far.We consider as the best setting the following: the primal heuristics, the four propagators, the OA cuts (for the epigraph model), and the standard branching rules that are part of SCIP, all enabled.The comparison is presented in Figure 10.As for the case k = 2, the dual bounds obtained with the epigraph model are better than the ones obtained with the quadratic model if plain SCIP is used.Although not many instances can be solved by using our techniques, the improvement that they bring to the solution process is significant in terms of primal-dual gaps.The instances solved to global optimality within the time limit are the smallest instances in our test set in terms of data points, i.e., instances 2 and 3 can be solved by the quadratic model and the epigraph model additionally solves instance 12; see Table 1 for details about these instances and Table 31 in the appendix for detailed results per instance.8.5.Discussion of the Numerical Results for Samples of Instances.The results discussed in the last two sections indicate that our techniques are highly beneficial for SCIP when solving the MSSC problem with a number of data points that is not too large.On the tested instances, this roughly means n < 1000 and d = 2.To further support this hypothesis, we conducted experiments on a broader test set of smaller instances.We have formed 10 new sub-test sets by extracting samples of sizes {100, 200, . . ., 1000} from the two-dimensional and large instances a1, a2, a3, s1, s2, s3, s4.Each sub-test set is comprised of 7 instances, which in total yields additional 70 instances.To sample the data, we use the Python routine skopt.sampler.Sobol; for Sobol sequences see Sobol' (1967).In Figure 11, we show two examples of these samples (or sub-instances).The blue points represent the original data points, while the red crosses represent the obtained sample.We believe that these samples extract meaningful information of the larger instances as the samples look similar to the full instances.In the following, we investigate how our novel methods scale with increasing problem size.We focus on the case k = 2, because the instances for k = 3 are still very challenging to solve and drawing reliable conclusions is difficult.
Primal Heuristics and Propagators.Table 12 shows aggregated results obtained by enabling our primal heuristics and propagators in SCIP.We can immediately see that SCIP's performance clearly improves if more of our components are enabled as more instances can be solved and the running times decrease drastically.In particular, enabling all our techniques allows us to solve using the quadratic model all 7 instances per sub-test set for up to 500 data points.Using the epigraph model, all instances with up to 900 data points can be solved to global optimality within the time limit.Moreover, for both models, we can solve almost all instances with 1000 data points if all our methods are enabled.Without our techniques, SCIP cannot solve a single instance even if the number of data points is 100.
Cutting Planes.We again only focus on the OA cuts and on the epigraph model.The summarized results are presented in Table 13.Interestingly, the OA cuts do more harm than good in this experiment.We can thus conclude that since these cuts are not based on problem-specific clustering ideas, they do not help as much as the propagators do for solving the MSSC problem effectively.
Branching Rules.Now we evaluate the performance of our branching rules.Again, we focus only on the two branching rules that yield some improvement in the solution process of the MSSC problem.The comparison results are displayed in Table 14.The entropy branching rule performs better if the epigraph model is used, whereas both the distance and the standard branching rules of SCIP are equally good if the quadratic model is used.
Best Setting.As in the last sections, we conclude the analysis by comparing plain SCIP with the best setting, see Figure 12.One can clearly see the significant improvement achieved by using our techniques.Moreover, it is also visible that as the sizes of the instances get larger, the harder the instances become.Therefore, if the size of the MSSC problem at hand is not too large, i.e., the number of data points n is around 1000, the dimension d is 2, and the number of clusters k is 2, then our techniques can be efficiently used to solve the problem to global optimality in just a few seconds.However, it is important to note that this is the case only for the instances we consider because the difficulty of a MSSC problem does not solely depend on the size of the instance but also depends on the structure of the given data points.
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.

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 "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-and-bound 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. (2021).Second, there further exist variants of the MSSC problem with additional side constraints as discussed in, e.g., Liberti and Manca (2021).Such side constraints allow for solution techniques that are feasibility-based, whereas all our techniques are optimality-based.Hence, a combination of both could yield an overall branch-and-bound framework that is even more effective for side-constrained MSSC problems.

Figure 2
Figure2.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 10 4 , then it assumes the gap limit of 10 4 in the plot.
Figure 3.Comparison of gaps using different propagators for 2 clusters.

Figure 4 .
Figure 4. Comparison of gaps using propagators with the OA cuts enabled or not for 2 clusters.

Figure 5 .
Figure 5.Running times and gaps comparison between plain SCIP and SCIP enabled with the best setting for 2 clusters.

Figure 6 .
Figure 6.Comparison of gaps using or not the heuristics for 3 clusters.

Figure 7 .
Figure 7.Comparison of gaps using different propagators for 3 clusters.
Figure 8.Comparison of gaps using different propagators and using or not the OA cuts for 3 clusters.
Figure 9.Comparison of gaps using different branching rules for 3 clusters.

Figure 10 .
Figure 10.Running times and gaps comparison between plain SCIP and SCIP enabled with the best setting for 3 clusters.

Figure 11 .
Figure 11.Example of samples extracted from the s1 instance.The sample sizes are: 200 (left) and 500 (right).The blue points represent the original data points, while the red crosses represent the sampled data points.

Figure 12 .
Figure 12.Running times and gaps comparison between plain SCIP and SCIP enabled with the best setting for the sampled instances and for 2 clusters.

Table 1 .
Information about the test sets.The first part corresponds to real-world test sets whose instances come from different sources.The second part corresponds to the synthetic test set.

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.

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.

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.

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.

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.

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.

Table 9 .
Comparison of mean number of nodes, mean running time per instance (in seconds), and number of solved instances using OA cuts for 3 clusters.

Table 10 .
Comparison of mean number of nodes, mean running time per instance (in seconds), and number of solved instances using or not the distance propagator for 3 clusters.

Table 11 .
Comparison of mean number of nodes, mean running time per instance (in seconds), and number of solved instances using different branching rules for 3 clusters.

Table 12 .
Comparison of mean number of nodes, mean running time per instance (in seconds), and number of solved instances using different propagators for the sampled instances and 2 clusters.

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.

Table 14 .
Comparison of mean number of nodes, mean running time per instance (in seconds), and number of solved instances using different branching rules.

Table 16 .
Results for 2 clusters using enabled heuristics.

Table 17 .
Results for 2 clusters using enabled heuristics and barycenter propagator.

Table 18 .
Results for 2 clusters using enabled heuristics and convexity propagator.

Table 19 .
Results for 2 clusters using enabled heuristics and convexity+cone propagator.

Table 22 .
Results for 2 clusters using enabled heuristics, propagators, cuts, and entropy branching rule.

Table 23 .
Results for 2 clusters using enabled heuristics, propagators, cuts, and distance branching rule.

Table 24 .
Results for 3 clusters using plain SCIP.

Table 25 .
Results for 3 clusters using enabled heuristics.

Table 26 .
Results for 3 clusters using enabled heuristics and barycenter propagator.

Table 27 .
Results for 3 clusters using enabled heuristics and convexity propagator.