Global Optimal K-Medoids Clustering of One Million Samples
Authors: Jiayang Ren
Kaixun Hua
Yankai Cao
Venue: NeurIPS 2022 Accept (2022)
Keywords: Large-Scale
Global Optimization
K-Medoids
Clustering
Lagrangian Relaxation
Branch and Bound
Bound Tightening

Abstract: We study the deterministic global optimization of the K-Medoids clustering problem. This work proposes a branch and bound (BB) scheme, in which a tailored Lagrangian relaxation method proposed in the 1970s is used to provide a lower bound at each BB node. The lower bounding method already guarantees the maximum gap at the root node. A closed-form solution to the lower bound can be derived analytically without explicitly solving any optimization problems, and its computation can be easily parallelized. Moreover, with this lower bounding method, finite convergence to the global optimal solution can be guaranteed by branching only on the regions of medoids. We also present several tailored bound tightening techniques to reduce the search space and computational cost. Extensive computational studies on 28 machine learning datasets demonstrate that our algorithm can provide a provable global optimal solution with an optimality gap of 0.1\% within 4 hours on datasets with up to one million samples. Besides, our algorithm can obtain better or equal objective values than the heuristic method. A theoretical proof of global convergence for our algorithm is also presented.


Global Optimal K-Medoids Clustering of One Million Samples
Jiayang Ren', Kaixun Hua Yankai Cao" University of British Columbia
Abstract
We study the deterministic global optimization of the K-Medoids clustering prob- lem. This work proposes a branch and bound (BB) scheme, in which a tailored Lagrangian relaxation method proposed in the 1970s is used to provide a lower bound at each BB node. The lower bounding method already guarantees the max- imum gap at the root node. A closed-form solution to the lower bound can be derived analytically without explicitly solving any optimization problems, and its computation can be easily parallelized. Moreover, with this lower bounding method, finite convergence to the global optimal solution can be guaranteed by branching only the regions of medoids. We also present several tailored bound tightening techniques to reduce the search space and computational cost. Exten- sive computational studies on 28 machine learning datasets demonstrate that algorithm can provide provable global optimal solution with an optimality gap of 0. 0.1% % within 4 hours datasets with up
million samples. Besides, algorithm can obtain better or equal objective values than the heuristic method. A theoretical proof of global convergence for our algorithm also presented.
Introduction
In this paper, we concentrate on the K-Medoids clustering problem: given a dataset, K-Medoids aims to select K samples from the dataset cluster medoids that minimize the of dissimilarities from other samples to the closet medoids. Here, selecting from existing samples as medoids is called the "medoids on samples" constraint. K-Medoids is similar the K-Means problem except for the arbitrary dissimilarity measurements and the "medoids on samples" constraint. Hence, K-Medoids generally has better interpretability than K-Means M.. K-Medoids problems are generally NP-hard to solve exactly. So there are many heuristic methods, such as PAM (2). CLARA [3], K-means-like method 4. and PAMM (5). However, none of them can guarantee reaching the global minimum.
To solve the K-Medoids problem deterministically. back in the 1970s, Coruuejols et al. 6 proposed a Lagrangian relaxation method for obtaining lower bound of the optimal value, which can guarantee a relative optimal gap smaller than 1/c. Interestingly, to the best of our knowledge, this is the only application in which the Lagrangian duality gap can be guaranteed within certain threshold. Several improvements to this method are proposed in the literature. One common direction is improving the heuristic process of Lagrangian dual problems. 7 constructed semi-I Lagi-Lagrangian relaxation, then accelerated the heuristic process by cutting plane method. It was further improved in [ by fixing variables in the heuristic process. Their works can deal with medium-scall datasets (c.g., 3000 samples). Another common direction is column generation 9
One limitation of these Lagrangian relaxation based methods is that there is guarantee of con- vergence to the global optimal solution, due to the existence of duality gap. To address the global
"corresponding author: yankai.car cataaiicaouulll @bc.ca. Jiayang Ren (rjy 12y1233076
and Kaixun Hua (kaixun hacixun.huaeee @ww.euuccaae contributed equally.
36th Conference on Neural Information Processing Systems (NcurIPS 2022).

---

optimum guarantee, few works have utilized the branch and bound (BB) scheme in K-Medoids [3. Christofides and Beasley [T]] proposed prototype by branching on binary variables and was capable of small-scale datasets (200 samples). This scheme was further improved by [12] through branching on semi-assignment constraints, which obtained solutions faster on small-scale datasets. [13] presented tighter form of the K-Medoids problem and was capable of medium-scale datasets (900 samples.. Besides these specific works, several off-the-shelf solvers also contain well-implemented BB schemes, such as CPLEX [44] and Gurobi [15. These solvers can typically address the K-Mediods problem with medium-scale datasets (e.g., 2000 samples).
The BB schemes used in the aforementioned works can only address medium-scale datasets because they need to branch on all the binary variables to guarantee convergence, which leads to poor scalability because the number of binary variables increases linearly as the number of samples. Therefore, we focus on the reduced-space BB scheme proposed in the stochastic programming community, which only needs to branch on the space of coupling variables (e.g., medoids). This idea was first reported for stochastic mixed-integer nonlinear programs in 177. The convergence for stochastic nonlinear programs was proved in [8]. Recently, the authors in [1]] adopted this method and reported the capacity of solving the K-means clustering problem (a similar clustering task without the "medoids on samples" constraint) with 210,000 samples in parallel to an optimal gap of 2.5%. However, the convergence of the reduced-space BB scheme can only be guaranteed when it is combined with suitable lower bounding methods. In previous works, lower bounds are constructed by the Lagrangian relaxation [16] 117 [9] removal [18] of the so-called aa-aiiccppatiiity constraints (an implicit constraint that all the samples share the same coupling variables, that is Constraint T0 in this paper). However, these methods unsuitable for the K-Medoids problem, since they cannot address the "medoids on samples" constraint (a constraint that has a large dimension and connects all samples) and thus cannot ensure the convergence of the reduced-space BB scheme. Moreover, these lower bounding methods require the solution of many expensive mixed-integer nonlinear programs to global optimality at each node. In contrast, the work of Cornuejols al. 6 generates lower bounds for the K-Medoids problem by the Lagrangian relaxation of the semi-assignment constraint (Constraint 2b in this paper), and the formulation directly takes the "medoids on samples" constraint into consideration. This method has a closed-fome solution, and its computation can be easily parallelized.
linear programs in [T6] and then extended to stochastic
Our Contributions: Firstly, we combine the method of reduced-space BB scheme and the work of Cornuejols et al. [6]. and adapt these methods for the K-Medoids problem. The Lagrangian relaxation in [6) is used to generate a tight lower bound, which already guarantees maximum gap at the root node. With this lower bound, we can prove the finite convergence (detailed proof' presented in the Appendix ( by branching only on the regions of medoids (the dimension of which is independent of the number of samples). Secondly, several tailored bound tightening techniques, including probing and feasibility based methods, are proposed significantly reduce the search space and computational cost of each node. Lastly, we propose an upper bounding method by heuristic updates of the candidate solutions obtained in the lower bounding procedure. These contributions enable our algorithm to be extremely scalable.
we provide open-source project in Julia and present numerical experiments on 28 datasets. Specifically, for two datasets with one million samples (which is 400 times larger than the state-of-art BB method), our algorithm can reach small gap (0.1%) 1) within one hour. To the best of our knowledge, no other article has ever reported results on datasets of this scale for K-Medoids problems with global convergence guarantee.
Capability For Million Scale Problems:
2 Preliminaries
K-Medoids problem: we denote a dataset with S samples and A attributes as X = #s}
SRSA in which x, = [I,,, aa R is the sth sample, and Fs,e is the ath attribute of sth sample. The K-Medoids problem aims to select a subset from the dataset as the cluster medoids with the following objective:
(.r,,.., "s.... (.s.....
min
Emin ||=, - 11 1 KEK.
WE M
#ES
where jut == lut, m E A is the medoid of kth cluster, μl := {u², ..., (² J the medoids set, S == {1, S} is the index set of samples in the dataset, K. := {1,.... ..., K) is the index set of clusters,
2

---

A:= (I, ...,. A} is the index set of attributes. We use E X to represent the "medoids on samples" constraint in the K-Medoids problem. To facilitate the discussion of the BB scheme, we introduce o M = {اسي - p S ji} as the initial regions of medoids. Here, ft and fi are the lower and upper bounds of medoids referred from data. Specifically, Pa = min 2,,0 and Fa = max a,a are the
sES
sE
bounds of the ath attribute. Note that the introduction of M does not affect the optimal solution.
Equivalently, the K-Medoids problem can also be expressed in an extensive form (EF) 20
(EF)
2(M)
min
ΣΣ
(2a)
b.y
SES jES
s.t.
Σ badd =1
(2b)
b..j S yj
(2d)
JES
bs.j yj E (0. 8, E S, k K
(2e)
Σ Uj K
(2c)
(2f)
JES
where d=.-=1|a,,- 1,113 the distance between samples I, and 2;; yj is equal to 1 sample z, is medoid, and otherwise O; bz, is equal to 1 if sample 2, belongs the cluster whose medoid is Zj.
and otherwise Note here ds.j is computed offline, so arbitrary dissimilarity measurements can be used without affecting the solution process. Constraint [ ensures that each sample v, must belong to one cluster. Constraint - guarantees that total of K samples are selected as medoids. Constraint 2 expresses that sample Vi be assigned sample æi only if 2g is the medoid of a cluster.
Lagrangian Relaxation of the K-Medoids problem: by the Lagrangian relaxation the Constraint pp. we can have relaxed problem of the Extensive Form 2
(M,
mine b.
X (day Aa)bas X,]}
(3)
JES
s.t. Constraints in the Extensive Form except
where A the Lagrangian multiplier for sth sample. We define the multiplier set A,, As}.
For a given A, (LD(M, A) can be solved analytically. Considering the objective function and constraints of the Relaxed Problem 3 note that the optimal values for ba.j given by
Yj if dajjjjj. Xa A..00 0,
(4)
10
otherwise.
We define the contribution of jth sample to the objective function as:
(A( Σ
min(0, daj -Aa).
(5)
SES
The Lagrangian problem can be reduced to
BLD(M
min 8 ()E(
Σps(A)
+ >
jES
#ES
(6)
Σ=5=K
Its closed-form solution be easily derived by selecting yj corresponding the K smallest Pj- Lagrangian dual: proper selection of may generate a tighter lower bound. Hence, we need to solve the Lagrangian dual problem get the tightest lower bound:
BLD (LD max (LD A).
There are several off-the-shelf heuristic methods update such as Sub-gradient method [T]] Volume method (1 and Cutting-Plane method (22). The solution to the Dual Problem provides
lower bound to the original problem, with (I.D < z(M). Particularly, for this problem, there is a unique guarantee maximum gaps between the lower bounds and the global optimal value, which
is z( (M) ( *, (
Σ max d=.j- This guarantee is theoretically proved in JES
23
3

---

we implement a parallel version of the tailored Lagrangian relaxation, in which the computations of contributions, Ø5 (A), are performed in each process. At the beginning of the BB procedure, given P processes, the index set S of the whole dataset is evenly divided into P subsets, Sp. p 1,.., P. Then, each process computes and stores the distance assigned it, which is d E S E Sp. At each BB node with M, denote the feasible index set of all the medoids as 3 Lagrangian Based Lower Bound
The aforementioned Lagrangian method cannot close the duality gap. To obtain global optimal solution, we propose to combine it with the spatial branch and bound (BB) scheme, in which Lagrangian relaxation provides a lower bound for each BB node.
Tailored Lagrangian Relaxation in BB nodes: Since our BB scheme branches on the regions of
the medoids μ4, at each BB node, we solve the problem with the partition set M C M. That is, constraint μ E M is added at each BB node. We further define == S 1 S n² medoid region of the kthh cluster, where 2º and j the lower and upper bounds of the ktth medoids. Since the possible region of each cluster's 's medoid can be different, we need to know the assignment of medoids to each cluster. Hence, we replace yjj with y in EF 2 vÿ is equal if sample zÿ is selected as the medoid of the kth cluster, and otherwise O. The Constraint 21 is replaced by
} as the
Σリソ K (2d) jES
ΣØ =1 (2d.1), jES
Σ REX
5 S (2d.2).
Notice that µ² is not directly modelled in EF, we need to transform the medoid regions μZ E Ma to the constraints of vf. The regions of medoids implicitly control the selection of y because have / 0, if ny & M Hence, we define the feasible index set of ktth cluster's medoid as (M) == [j E S 1 < xx < j² In this way, we can reformulate the reduced Lagrangian Problem at the BB node with M as follow:
BLD(M, A)
min JEE (M)
P5(A) + >EE Σ
REK. (ES++(M)
sES
(8)
s.t. Constraint (2d.1), (2d.2).
with the Lagrangian dual problem defined as BLD (M) max BLD (M, A). The optimal value of y (2) let
in Problem & can be obtained by the following procedure: (1) initialize h and
jf be the index of smallest Pم (A) subject to S (M) and J"; (3) set 2/ 1 if j otherwise, 0; (4) add jK* J"; (5) update k=k + and return to (2) until K.
Parallel implementation:
S+(M)= = S (M) US U.. US* (M).
Then, each process computes the assigned contributions using the stored ds.j. which is (A) Vj E Sp (MM. The remaining parts of the parallel implementation (e.g., update of X) are identical to the serial version.
4 Tailored Reduced-space Branch and Bound Scheme
This section introduces the tailored seduced-space branch and bound algorithm for the K-Medoids
problem, which is detailed Algorithm Specifically, starting from the root node with M, the BB scheme recursively partitions the medoid regions. At each iteration, the algorithm first selects the node with the lowest lower bound and denotes it as M. Then, bound tightening techniques are applied to get a tightened medoid region M. The lower and upper bounds of the tightened node are computed and updated sequentially. the gap larger than the tolerance, we branch on the medoid variable p with the maximum range to get two sub-nodes Mj and M2 with relint( M2) n relinf( = 0. where relint (M) the relative interior of M. On the convergence of the algorithm, we have the following theorem with detailed proofs in Appendix A
Theorem Algorithm|y / convergent to the global optimal solution after a finite step L,, with 3r a,, by only branching on the space of p.

---

Algorithm 1 Tailored Reduced-space Branch and
Algorithm 2 Probing
Bound Scheme for K-Medoids clustering
Select region the branched variable M
1: Initialization
2: Divide M into two sub-regions Msub.
2: Load and store dataset as X.
3: Compute basic LB for cach sub-region.
3: Initialize the iteration index 1 = 0, the node set
4: Delete the sub-region if Brasse (Maub) > ou-
M M, and tolerance 2 0.
Algorithm 3 Cluster Assignment
4: Compute initial lower, upper bounds: 31 3(A)).
for Sample æ. n M, do if wt is not assigned in the parent node then if d' d' < da Vk E AC) ( k* then
on = =.MM..
2:
5:
while Mi & and | anl 2 do
3:
6:
Node Selection
Assign v, to k* cluster.
7:
Select and delete from Node M E M with its parent node's lower bound [ = & (for the root node, select and delete itself).
5:
end if
6: end if
7:
end for
8:
Update / 1
Bound Tightening
Algorithm Feasibility Based Bound Tightening
10:
Probing Alg.
1: for Cluster k E K. do
11:
Cluster Assignment // Alg. 3
2: Obtain the assigned set Shi according to Alg
12:
Feasibility Based Bound Tightening // Alg. 4
3: Compute RHS cx -B(M)
13:
Tighten according the "medoids samples"
for Attribute E do 5: Obtain the tightened node, M - (n
constraint.
14:
Obtain the tightened node M.
1/a S f.o }. by solving the incquation
15:
Bounding
16:
Compute lower, upper bounds: (MM)) (M))
6:
Update M M n M
7: end for
17:
Update & + (nin(((M)) M' E MI.
8: end for
18:
Update dy + min min(ai-1.0 a (M)]-
19:
Delete all MM' from M if ((M)) > au.
Algorithm 5 Lagrangian Based Lower Bounding
20:
Branching
1: Compute basic LB, (M).
21:
[M n X] > and 18t -0i-oal 0z 2 e then
2: if [Sbasic (M) axl 2 e then
22:
Partition M into subsets M, and M2 with
3: Compute Lagrange Dual LB, (LD(M).
relint (M))
) n (elint(M2) =
23:
Add M, M. if M" n. + 0, Vk E i
4: B(M) mar(Bbasie (LD).
{1,2}.
5: else
6: ((M) (M))se(M).
24:
end if
25:
end while
7: end if
4.1 Lower Bounds
Basic lower bound: Besides the Lagrangian based lower bound, we also introduce a basic lower bound Bousic- Although it is not as tight as the Lagrangian based lower bound BLD. its computational cost is significantly lower. Thus it is useful when the BB tree is deep with small medoid regions. For K-Medoids I there an implicit constraint that all the samples share the same medoid set:
1. - Vs (1,..., 1), E K.
(10)
which uē is the kith medoid for the sth sample. By removing the "medoids on samples" constraint E X and Constraint 10 we have the basic lower bound (hasic(M) with Biusio (M) S =(M):
Brasio (oeico(M (M)
min Σ minn ||,, 1 112 = Σ min min 1|,, -1, 11ll.. kER
(11)
P.EM
KEE MEEM*
SES
SES
We define df mi (M) == min ||2, - 1 1 and obviously Prasic(M)
Σ min di SE KEK mmmin (M). The
analytic solution to dls mim (M) can be easily derived: / mid Fa,as at), a E A, in which _h and 3 are the bounds the ath attribute in kth medoid region.
Lagrangian based lower bound: as shown in Algorithm 5 by combining the basic and LD lower
We note that because (M) S =(M) and (r.n(M) < (MM) we have ((M) S x Besides, according to the guarantee of maximum gap between the Lagrangian based lower bound and the global optimal value, we can also derive guarantee of maximum gaps for ((M):
bound methods, we give the final form of lower bounds: 3(M) mar{Bbasic (LD
Proposition 2. =(M)-8(MM
e. where 2,= Σ max ds.j-
-x-B(M)
SES JES
5

---

4.2 Upper Bounds
Upper bounds of the K-Medoids problem can obtained by selecting a feasible set of medoids M n Given feasible solution, we can obtain the assignments of samples and object value through Equation I without solving any optimization problem explicitly. However, the selection jr will critically influence the efficiency of convergence. Hence, in implementation, introduce two efficient upper bounding methods:
we use the K-Means-like method 4 as backbone. Because the K-Mcans-like method sensitive to the initial seeds and apt to the local minimum [5). we run this method several times with random seeds to get a candidate solution set. Then, we improve the quality of the upper bound by feeding the candidate solution set as initial guesses to Evolutionary Centers Algorithm [24].
Root node:
Child nodes: we utilize the candidate solution obtained in the lower bounding method to update the upper bounds. Specifically, we first compute the corresponding assignments, then feed these assignments as the initial solution to the K-Means-like method. It should be noted that this K- K-Means- like method gives solution full-filling aill-filling the "medoids on samples" constraint in K-Medoids.
Bound Tightening Techniques
Here, we introduce some bound tightening techniques tailored for the K-Medoids problem to reduce the search regions M at each BB node while guaranteeing the optimal solution is not excluded. These bound tightening techniques are deployed before lower and upper bounds calculation.
5.1 Probing
Probing is a technique reduce bound regions by exploiting the inequalities in mixed integer programming problem 251 This paper exploits the implicit inequality that the lower bound of a region containing the optimal solution is smaller than the current best upper bound. As shown in Algorithm 2 it tentatively restricts the bounds of μ. to subinterval and then computes the corresponding lower bound. If the lower bound is larger than the best upper bound, can conclude that no optimal solution exists within the subinterval, and tighten the bound on the variable accordingly. The computation costs are trivial since we use the basic lower bounding method.
5.2 Feasibility Based Bound Tightening
This subsection introduces the Feasibility Based Bound Tightening (FBBT) technique. It is based the observation that the assignment of some samples to clusters can be predetermined.
M2
dlastn
damis
Figure 1: Cluster assignment with 3 clusters. In this example, we have df dzymaz dz
d)
s,min and
d3 mim- Therefore, we assign Fg to the second cluster.
Cluster assignment: if we have ||r, n** | ||x, n²* 122 Vk- E ( then obviously sample vs is in the 1+thh cluster. However, the value of μ here is not known before solving the
overall problem. Recall that have mi (M):= min ||=r, ~ Similarly, also define
=*E
dh Msm (M):= max ||x, 11 112 the maximum distance between F, and μa with μ M². At WM
a BB node with partition M, if we have de'mna (M) < dh demmin (M), Vk- EK),t,, then sample check whether
guaranteed to be assigned to the k th cluster. Based this observation, the cluster of some samples can be determined a BB node as shown in Algorithm| 3 Recall that
d de.min (M) is already computed as byproduct of the basic lower bound in Section Similarly, since the cluster regions for each node are rectangles, d (M) can also analytically computed. Besides, since child nodes are subsets of the parent node, once sample is assigned to cluster in

---

parent node, it is also assigned to the same cluster in its child nodes. Figure shows a simple
illustration of this cluster assignment scheme.
FBBT: In a node with M, we denote the index sets of samples assigned to cluster k as S FBBT infers the regions of 2 from the following inequation:
3(M)+ +
Σ - =
(12)
SESA
where a is the current best upper bound, ( represents the lower bound contributions of samples unassigned cluster k (either samples assigned other clusters or samples whose assignments are undetermined), Similar to the calculation of the basic lower bounding method, we have (MM)
Σ min d (..mmn SES)SS RER
As shown in Algorithm 4 FBBT applied for each cluster. Specifically, for the kth medoid region. 1) obtain S according to the results of Algorithm 3 2) compute the fight-hand-sid RHS
ag ((M)=; 3) for each attribute E A. obtain the tightened region, = - p P
by solving the following inequation:
Σ (Ix.a Ma) 2 + Σ Σ
min (#.. - 11)² RHS, EeE E. SESA :EAA "TEM
(13)
RESA
Since 13 a simple quadratic inequation, we can analytically solve the problem to obtain the lower and upper bound of f as L. and f.a Finally, update M by M + M n M
5.3 Miscellaneous
Symmetric-breaking: we enforce a cammetric--reakinng constraint on the first attribute, i.e., 1 S "1 k+1 Vk = I,..., K 1. It is used to tighten the regions of M at each BB node.
Bound tightening according to the "medoids on samples" constraint: since the K-Medoids problem requires medoids be selected from samples the dataset, we can further tighten the regions of M by M- {pl A S μ S (). *EM
min Zs,as and jtc
max Fs.a. .EMпX
5.4 Effects on Computation
With reduced search space, these bound tightening techniques can significantly reduce the number of BB nodes to explored. Moreover, they can also accelerate the calculation of bounds. For example, for the Lagrangian Problem 8 we need to calculate the contributions, Pj. of samples belonging the overall feasible set of medoids, S+ $1+ U 52+ U... SK++ A reduced search space MM means reduced overall feasible selection set, which further leads to fewer calculations and faster speeds.
Experiments and Discussions
To evaluate the performance, we use 25 datasets from the UCI Machine Learning Repository [26). one dataset called PR2392 from (7 one dataset called HEMI and one synthetic dataset from
in the clustering category with multivariate numerical attributes and range from 100 to 2,458,285 instances. We implement our branch and bound K- Medoids algorithm with Lagrangian based lower bound method (BB+LD) in Julia 1.6. 1. We also implement four other methods to compare the performance, including the extensive form of K- Medoids problem solved by the state-of-art global optimizer CPLEX 20.1.0 14] (CPLEX), the heuristic method in the root node (Heuristic), the pure Lagrangian relaxation method without branch and bound (LD), and the of BB method for K-Means in [19 adapting to K-Medoids by adding "medoids on samples" constraint (BB+Basic). We executed all the experiments on high- performance computing cluster, of which each node contains 40 Intel cores at 2.4 GHz and 202 GB RAM. We present the results of cluster number K = 3 in the main body with additional The complete code files can be found https
19 The datasets from UCI Repository
results K = and K 10 in Appendix C
//gi frgtthub.crm com/ wom/YankaiGrroup//ooooooooaoaaaa wwtpoouuugggoollllllmomoooo kmedoids

---

The numerical performances are compared using three indicators, including upper bound value (UB), relative optimal gap, and the number of solved nodes. Here, UB the best upper bound in all the iterations. The relative optimal gap is the relative difference between the best lower and upper bounds, which is calculated by Gap(%)= = ou on-0t
min(ar.8) x 100%. This relative gap represents the worst gap of the found solution from the global optimal solution. Notably, this gap is a particular property of deterministic global optimization algorithms. Heuristic methods can not quantify the worst gap from the global optimal solution. The number of solved nodes represents the total number of BB nodes.
All the deterministic methods CPLEX, LD, BB+Basics, and BB+LD, share the same terminal criteria, including (1) Gap < 0.1% (2) the solving time reaches 4 hours, (3) the number of solved nodes reaches million. Both LD and BB+LD use the same sub- robdgradient method called the Volume method [21. Particularly, for BB+LD, the maximum iteration of Lagrangian multiplier updates in each node 10. For LD, the maximum iteration is 10,000. In the heuristic method, we execute the K-Means- like method 20 times with fixed seed set including 20 different seeds, then feed the solution set as initial guesses to ECA to get the best solutions.
Table 1: Serial numerical results (K
DATA- SET
SAM-
PLE
Dis- ENSIGN
METHOD
UB
NODES
GAP (%)
TIME: (S)
DATA- SET
SAM-
PLE
DIM- ENSION
METHOD
UB
NODES
GAP
(%)
TIME:
(8)
HRURISTIC
HEURISTIC
9.91E+06
CPLEX
83.91
20.10
41
CPLEX
9.91E+06
$0.10
2044
IRIS
150
LD
83,91
1
20.10
19
НEMI
1955
LD
9.91E+06
12.45
3176
BB-BASIC
83.91
2.0E+5
50.10
124
BB+BASIC
9.91E+06
1.9E+6
4.32
4H
BB+LD
25
20.10
BB+LD
20.10
97
DEURISTIC
$98.29
HEBRISTIC
2.13E+10
CPLEX
CPLEX
2.13E+10
$0.10
5339
SEEDS
210
2.13E+10
35
BB+BASIC
598:29
2.4845
20.10
273
BB4BASIC
2.13E
3.0E45
4.38
411
BB-LD
598.29
9
$0.10
84
BB+LD
2.13R-10
37
$0.10
123
HEURISTIC
HEURISTIC CPLEX
CPLEX
629.02
20.10
60
GLASS 214
LD
629.02
1
20.10
19
TRR
5456
24
LD
1995
0.41
4n
RR-RASIC
2,2946
2871
BRABASIC
1996+005
4F44
1644.88 88
BB+LD
629.02
32
20.10
107
BE+LD 1.9E+05
553
20.10
325
[EERRISIC
8.65E+05
HEURISTIC 2211.19
CPLEX
8.632+05
CPLEX 2199.12
249
20.10
7195
50.10
287
BB-BASIC
8.63E+05
5.0E+6
25.74
7698
BB+BASIC BB+LD
2200.04
1. 1.34+5
573.92
4H
BB-LD 8.63E+05
53
86
2199.10
50.10
222
HEURISTIC
CPLEX
50.77
50.10
74
CPLEX
UK
258
LD
50.77
0.96
176
RDS_
10000
LD
1,49+007
14,89
4H
9.5
2.68+5 1,49+007
5883
BB-LD
50.77
119
20.10
89
BE+LD
31
20.10
203
HIURISTIC 7. 2.83E+11
HOURISTIC
8.21
7833+11
299
7
17898
82114077
22,97
411
BB-BASIC
7.83-+11
9.6E+4
$0.10
173
BB+BASIC
8.21E+07
2.36+5
25.48
4H
BBeLD
7.83E+11
99
107
BB+LD
8.21E+07
465
$0.10
1555
HRURISTIC
34
HEURISTIC 1.95+007
CPLEX
8.33--10
50.10
124
CPLEX LD 1.95+007
WHO
440
LD
1.33-+10
11.10
37
36733
1.43
4H
8.33E+10
5.0E+6
11943
1.95++07
4.3E+4
298.13
4k
BBaLD
8.538+10
193
20.30
117
BB+LD
1995+007
101
20.10
1936
HIURISTIC
2 25755406
HOURISTIC 486.75
2.75E+06
CPLEX LD
HCV
572
LD
2.75E+06
0.34
536
RDS
50000
486.75
4.34
4H
BB+BASIC
2.75E:06
1.9E16
20.10
5587
BB4BASIC
476.79
2.0845
35.72
4t1
2.75F406
4818
811
HEURISTIC
2.625406
HEURISTIC CPLEX
4,94F+08
CPLEX
2.62E+06
23
0.23
285
ABs
740
2.62E+06
726
KEGG
53413
4.94E+08
151.02
4nt
BR-RASIC BB+LD
2.62F406
2.4F46
57.45
4,044008
1.8444
11.72
411
2.62F.:06
125
20.10
119
BR+LD
4.94R-08
177
50.10
3901
DEERISTIC 1136.93
HEURISTIC 1,26E+05
855
TR
980
BB+BASIC
L.D
1136.93
1136.93
2,1646
209 23
0.80
130
4u
URBAN
GB_10
100000
BB4BASIC
LD
1.15E+05
1,26E+05
7.68:4
17 17.86802
26.5704
4n1
4n
1.15K-05
6834
HEURISTIC
1,28E+09
CPLEX
1
$0.10
347
SGC
1000
1.28E:09
21.85
1.28F409
3.8F45
2739
BB+LD
1.28E+09
189
20.10
140
"OUT-O--MEMORYY
8

---

6.1 Serial Results
Table [ shows the serial numerical results for the datasets with less than 100,000 samples. In terms of best found upper bound (UB), the deterministic global optimization methods (CPLEX, BB+Basic, and BB+LD) can obtain better UB than the heuristic method in several datasets (e.g.. IRIS, BM, RDS). While try to obtain the best heuristic results by multiple initializations, there is still no guarantee for the heuristic method to reach the global optimum. Moreover, the heuristic method can not provide an optimal gap estimate the solution quality. Notably, method and CPLEX provide the same upper bounds and gaps for the datasets smaller than PR2392 (2392 samples).
As for the gaps and running times, CPLEX and LD can obtain a small gap (50.1%) with less time when datasets are smaller than HCV (572 samples). However, when datasets are larger than HCV, our method, BB+LD, can reach the small gap with much less time than CPLEX and LD. Remarkably, our method can reach the small gap all the datasets up to 100,000 samples within 4 hours. CPLEX fails to give a solution because of the out-of-memory error in the root node when datasets are larger than PR2392 (2392 samples). L.D and BB+Basic methods may not reach the small gap within 4 hours or maximum iteration limits when datasets are larger than BM (249 samples). Recall that LD has no guarantee of reaching the global optimal solution because of the duality gap. All these results show that our method is extremely scalable for K-Medoids problems even under serial mode.
Notably, we followed the tuning guide in the CPLEX documentation improve its performance. For example, we initiated the root node with the Heuristic solution and specified the node file and memory parameters. Moreover, we tried different MIP strategies, including start, pre-solve, and heuristic strategies. In T we present the best results obtained by CPLEX.
6.2 Parallel Results
For BB+LD, we executed parallel experiments utilizing MPI interfaces in the large and huge scale datasets with K = 3. The results of large-scale datasets ranging from 10,000 100,000 samples presented in Appendix B Table Sssows the parallel numerical results for huge datasets with 100,000 to 284 samples. In this table, even for the dataset with 1,046,910 samples, BB+LD can still converge to the global optimal solution (50.10 10 gaps) within one hour, which is 400 times larger than the size of problem that the state-of-art global solver CPLEX can address. For USC1990 with 2,458,285 samples and 68 features, BB+LD can obtain a 6.33% gap within 4 hours. According to Theorem we can obtain global cooptimal solution for USC1990 given more calculation time.
Table 2: Parallel results of huge scale datasets (BB+LD, K 3)
DATASET
SAMPLE DIMENSION
CORES
UB
NODES
GAP(S) TIME(S)
RNG_AGR
199,843
7
1600
8.23E+11
99
$0.10
341.0
URBANGB
360.177
2
1600
4.14E+05
57
<0.10
327.0
SPNET3D
434.874
3
1600
2.28F++07
115
50.10
865.0
RETAIL
$41,909
1600
6.80E+09
20.10
80.0
SYNTHETIC
1,000,000
6000
9.44E+06
3
$0.10
171.0
RETAIL-II
1,046.910
2
6000
2.31F+10
214
50.10
2515.0
USC1990*
2,458,285
68
3000
6.91E+08
25
6.33
4H
daj COMPUTED ON-THE-FLY, NOT PRECOMPUTED STORED.
7 Complexity Analysis
7.1 Time Complexity
Denoting the number of samples S, attributes A. and clusters as K, our algorithm branches on the medoid regions to guarantee convergence, which contains A x K variables and is independent of the number of samples S. However, the number of BB nodes to converge is hard to predict since the efficiency of bounding methods depends on datasets. Also, the number of Lagrangian Dual (LD) iterations in each BB node is different because we set dynamic stopping criteria in case there is update of LB in several continuous LD iterations. The time complexity of one LD iteration is ((52) when d are precomputed. If dæ.j are not precomputed, the time complexity is O(A8)).
Table compares computing time and nodes on synthetic datasets. In this table, the average run- time of one LD iteration remains almost the same as the dimension changes and increases almost

---

quadratically as the number of samples increases. As for the BB nodes, although the number of branching variables is A x K. the actual number of explored BB nodes does strictly follow the growth of dimensions. These results are consistent with the complexity analysis mentioned above.
Table 3: Comparison of time and nodes on synthetic datasets with varies dimensions and samples
SAMPLES DIMENSION
20
30
40
50
60
70
100
NODES
10
3
3
31
19
3
3
39
69
1,000
TIME (MS)
4.28
4.28
4.06
3.57
1.62
4.02
4.11
3.78
4.32
3.8
3.21
3.79
17
15
10,000
TIME (MS)
423
363
387
357
390
361
373
366
369
363
211
358
100,000
NODES
7
27
31
3
15
27
13
25
21
35.2
28.7
35.6
37.1
NUMBER SOLVED NODES TIME CONSUMED IN EACH LAGRANGIAN DUAL ITERATION.
PROCEDURE CONVERGENCE.
7.2 Space Complexity of Distance
When distances dæ are precomputed and stored, the space complexity of ds.j matrix is O(5²). In our implementation, the datasets with no more than 100,000 samples were computed on one compute node with 40 cores and 202GB RAM. For more than 100,000 samples, we executed the experiments on multiple compute nodes. Each core precomputes and stores part of the da.j matrix, as described in Section 3 For example, we executed the one-million dataset RETAIL-II on 150 nodes with 6,000 cores. Each core only needs to store 1.25GB of da, matrix. For more than 1,000,000 samples, we calculated dæ.j on the fly, without precomputing and storing. this case, the time complexity of one LD iteration increases from 0(5²) to O(42²). Hence, we expect an acceptable slowdown when the dimension of datasets is small. For example, in Table the dataset USC1990 with two-million samples and 68 features was executed without precomputing dz.j+ and the result is acceptable. Table 4 compares the computing time of precomputed distance and on-the-fly computing. We can conclude that although calculating ds.j on the fly is slower than precomputing da.j. the slowdown is acceptable.
Table 4: Comparison precomputed distance and on the fly computing
TOTAL
TIME PER
DATASET DIMENSION
RUN TIME (8)
LD ITERATION (MS)
ON-THE-FLY PRECOMPUTED ON-THE-FLY PRECOMPUTED
ABS
21
34
12
1.42
0.70
HEMI
7
17
13
9.91
2.02
RDS RD8_ CCT
47
30
175
37.80
TR TRR
2.54
0.73
24
110
10.80
8 Conclusion
We presented a scalable global optimization algorithm of the K-Medoids problem by applying a tailored reduced-space spatial branch and bound scheme. This algorithm includes Lagrangian based lower bounding method and basic lower bounding method. Bound tightening techniques are also proposed to accelerate the solution process. We demonstrate algorithm's 's scalability by extensive numerical experiments and prove the convergence by theoretical analysis. Besides, one interesting finding that our algorithm can reduce the objective values for some datasets while remaining the same as the heuristic method for most datasets. We also compared other clustering evaluation metrics in Appendix D such Normalized Mutual Information (NMI) and Adjusted Rand Index (ARI). A similar conclusion with objective values holds for these metrics. However, it should be noted that only with the global optimal solution obtained in our work can we give a fair comparison of the heuristic methods and, in turn, contribute to developing a more efficient heuristic method.
Acknowledgments and Disclosure of Funding
The authors acknowledge funding from the discovery program of the Natural Science and Engineering Research Council of Canada under grant RGPIN-2019-0540 and the computing resources provided by SciNet (www.scinetie www.sccnnthppcca. and Digital Research Alliance of Canada (wwwaalliancecan.c)) Jiayang Ren acknowledges the financial support from the China Scholarship Council.
10

---

References
[1] Adil Fahad, Najlaa Alshatri, Zahir Tari, Abdullah Alamri, Ibrahim Khalil, Albert Y. Zomaya, Scbti Foufou, and Abdelaziz Bouras. A Survey of Clustering Algorithms for Big Data: Taxonomy and Empirical Analysis. IEEE Transactions on Emerging Topics Computing, 2(3):267-279, 2014.
[2] Leonard Kaufman and Peter J. Roussceuw. Partitioning Around Medoids (Program PAM). In Finding Groups in Data, pages 68-125. John Wiley & Sons, Ltd, 1990.
[3] Leonard Kaufman and Peter J. Rousseeuw. Clustering Large Applications (Program CLARA). In Finding Groups in Data, pages 126 163. John Wiley & Sons, Ltd, 1990.
[4] Hae-Sang Park and Chi-Hyyck Jun. A simple and fast algorithm for K-medoids clustering. Expert Systems with Applications, 36(2):3336 (13333333333,,,,,,, 2009.
[5] Erich Schubert and Peter J. Roussceuw. Fast and cager k-medoids clustering: O(k) runtime improvement of the PAM. CLARA, and CLARANS algorithms. Information Systems, 101:101 804, 2021.
[6] Gerard Comuejols, Marshall I.. Fisher, and George I.. Nemhauser. Location of Bank Accounts to Optimize Float: An Analytic Study Exact and Approximate Algorithms. Management Science, 2388):88--800, 1977.
[7] C. Beltran, C. Tadonki, and Ph. Vial. Solving the p-Median Problem with a Semi Semi-Lagrangian Relaxation. Computational Optimization and Applications, 35(2):239-660, 2006.
[8] Andréa Cynthia Santos. Solving large p-median problems using Lagrangean heuristic. Technical report, Université Blaise Pascal (Clermont Ferrand 2). 2009.
[9] Sergio García. Martine Labbé, and Alfredo Marín. Solving Large p-Median Problems with Radius Formulation. INFORMS Journal on Computing, 2345)5546-556, November 2011.
[10] Anton V. Ushakov and Igor Vasilyev. Near-optimal large-scale k-medoids clustering. Information Sciences, 145:344-362, 2021.
[11] N. Christofides and J.E, Beasley, tree search algorithm for the p-median problem. European Journal of Operational Rexearch, 10(2): 196-204, 1982.
[12] Alberto Ceselli and Giovanni Righini. A branch-and- anccc-eed--riceeeeee algorithm for the capacitated p-median problem. Networks, 45(3): (253))122-112 142, 2005.
[13] Sourour Elloumi. A tighter formulation of the p-median problem. Journal of Combinatorial Optimization. 19(1):69-88 83,, 2010.
[14] IBM IL.O Cplex. V20.1.0: User's Manual for CPLEX. International Business Machines Corporation, 2020.
[15] Gurobi Optimization, LLC. Gurobi Optimizer Reference Manual, 2022.
[16] Claus C Carée and Rüdiger Schultz. Dual decomposition in stochastic integer programming. Operations Research Letters, 24(1-2):37 1999.
[17] Ramkumar Karuppiah and Ignacio E Grossmann. A lagrangean based branch-and-ca autt algorithm for global optimization nonconvex mixed-integer nonlinear programs with decomposable structures. Journal of global optimization, 41(2):163-186 2008
[18] Yankai Cao and Victor M. Zavala. A scalable global optimization algorithm for stochastic nonlinear programs. Journal of Global Optimization, 75(2):393 416, 2019.
[19] Kaixun Hua, Mingfei Shi, and Yankai Cao. A Scalable Deterministic Global Optimization Algorithm for Clustering Problems. In International Conference Machine Learning, pages 4391-4401. PMLR, 2021.
[20] Hrishikesh D Vinod. Integer programming and the theory of grouping. Journal of the American Statistical association, 64(326):50653 519, 1969.
[21] Laureano F. Escudero, M. Araceli Garín, Gloria Pérez, and Aitziber Unzueta. Scenario Cluster Decompo- sition of the I Lagrangian dual in two-stage stochastic mixed 0-1 optimization. Computers & Operations Research, 40(1) (661)736 2013. 00049.
[22] Kibaek Kim and Victor M. Zavala. Algorithmic innovations and software for the dual decomposition method applied to stochastic -ixed-inteegr programs. Mathematical Programming Computation, 10(2):225- 266, 2018

---

[23] Marshall L. Fisher. The Lagrangian Relaxation Method for Solving Integer Programming Problems. Management Science, 50:1861-1871, 2004.
[24] Jesús-Adolfo Mejía-de Dios and Efrén Mezura-Montes. New Evolutionary Optimization Method Based Center of Mass. In Kusum Deep. Madhu Jain, and Said Salhi, editors, Decision Science in Action:
Theory and Applications of Modern Decision Analytic Optimisation, Asset Analytics, pages 65-74. Springer, Singapore, 2019.
[25] Savelsbergh. Preprocessing Probing Techniques for Mixed Integer Programming Problems. ORSA Journal on Computing, (4):45 454, November 1994.
[26] Dheeru and Casey Graff. UCI machine learning repository. 2017.
[27] Manfred Padberg and Giovanni Rinaldi. branch-and-cut algorithm for the resolution of large-scale symmetric traveling salesman problems. SIAM review, 33(1):60-1 100, 1991.
[28] Edward Wang, Genpei Cai, Riley Ballachay, Yankai Cao, and Heather L. Trajano. Predicting Xylose Yield Prehydrolysis of Hardwoods: Machine Learning Approach. Frontiers Chemical Engineering. page 84, 2022.
[29] Reiner Horst and Hoang Tuy. Global optimization: Deterministic approaches. Springer Science Business Media, 2013.
12

---

Checklist
1. For all authors...
(a) Do the main claims made in the abstract and introduction accurately reflect the paper's
contributions and scope? [Yes]
(b) Did you describe the limitations of your work? [Yes] See Section ??
(c) Did you discuss any potential negative societal impacts of your work? [N/A] Our work is about a global solver of the K-Medoids problem. The problem itself may have societal impacts.
(d) Have you read the ethics review guidelines and ensured that your paper conforms to them? [Yes]
2. If you are including theoretical results...
(a) Did you state the full of assumptions all theoretical results? [Yes]
(b) Did you include complete proofs of all theoretical results? [Yes] See Appendix
3. If you ran experiments...
(a) Did you include the code, data, and instructions needed reproduce the main ex- perimental results (either in the supplemental material or as a URL)? [Yes] See the supplemental material
(b) Did you specify all the training details (c.g.. ... data splits, hyperparameters, how they were chosen)? [Yes] See Section and the supplemental material
(c) Did you report error bars (e.g.. with respect to the random seed after running experi- ments multiple times)? [N/A] Our method is deterministic method.
(d) Did you include the total amount of compute and the type of resources used (e.g., type of GPUs, internal cluster, cloud provider)? [Yes] See Section 6
4. If you using existing assets (e.g., code, data, models) or curating/releasing
assets...
(a) If your work uses existing assets, did you cite the creators? [Yes] We compare our
method with some existing methods and we cite the creators in the main body.
(b) Did you mention the license of the assets? [N/A]
(c) Did you include any new assets either the supplemental material or as URL? [Yes] We include our project codes in the supplemental material.
(d) Did you discuss whether and how consent was obtained from people whose data you'l 're using/curating?? [N/A] We use the datasets from UCI Machine Learning Repository and three open-source datasets.
(e) Did you discuss whether the data you are using/curating contains personally identifiable information offensive content? [N/A] The meaning of datasets are described in UCI Machine Learning Repository.
5. If you used crowdsoureing or conducted research with human subjects... [No] Our work does not use crowdsourcing or conducted research with human subjects
(a) Did you include the full text of instructions given to participants and screenshots, applicable? [N/A]
(b) Did you describe any potential participant risks, with links Institutional Review Board (IRB) approvals, if applicable? [N/A]
(c) Did you include the estimated hourly wage paid to participants and the total amount spent on participant compensation? [N/A]
13


Reviews:
[{"review_id": "sN7nBna7Mmf", "reviewer": "Reviewer_DaVW", "summary": "\nThe authors propose a branch and bound approach to solving the K-mediod problem by taking advantage of the Lagrangian relaxation method explored by Cornuejols in\n the 70's. In particular, they modify the lagrangian relaxation for each cluster by building appropriate boxes for the BB scheme and provide various speedups by\nprobing, feasibility based bound tightening and parallelization. Experiments are shown on a wide variety of datasets.\n", "questions": "(A) How does the scaling look like in terms of the dimensions d? Also, in general what is the computational complexity of the proposed algorithm? Of course, it w\nould depend on the number of the branch and bound steps. \n\n(B) The algorithm seem to be more sensitive to dimensions than to samples in terms of performance (Table 2). Is that really true and if so, any reason why that could be the case?\n", "limitations": "The key reason it works is because of the Lagrangian relaxation available to the K-mediod setting and it is not clear if this will extend to other clustering problems. Can the authors shed some light as to where else this could be applicable?", "rating": 7, "confidence": 4, "soundness": 4, "presentation": 4, "contribution": 3, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "The paper is clearly written with references to related work and is easy to follow. The theoretical results are shown in the main paper with the proofs delegated\n to the supplementary. Additional experiments with varying sizes of clusters are also provided.\n\nThe problem is of broad interest since clustering is a core component of many machine learning tasks. The key ideas of using Lagrangian relaxtion and branch and\nbound techniques have been explored in the broader literature but combining them in the context K-mediods problem along with mediods on samples constraint makes\nit novel and relevant.      \n\nThe paper pretty much builds on the work of Cornujols et al as well as the branch/bound literature and the proofs of finite termination/global optimum are provided.  Most of the literature for the relevant work is covered and referenced in the paper.\n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "FqJtFFnhIB", "reviewer": "Reviewer_2o3M", "summary": "The paper presents a spatial branch-and-bound method for the K-medoids clustering problem. The authors introduce several bound tightening techniques tailored for their approach. In numerical experiments they demonstrate that their method is able to solve large scale instances of the problem to global optimality.", "questions": "The presentation of the results is mostly clear, but I have a few questions/comments:\n\nIn lines 96-97 you are mentioning that the value $d_{s,j}$ are precomputed. However, the number of values scales quadratically with the size of  the data set, which seems intractable for the claimed size of 1M samples. How do you deal with this issue?\n\nIn line 130, the authors define $S^{k+}$ as the set of indexes whose samples are contained in the $k$-th medoid constraint set. Then, in equation (9) we consider the union $S^+ = S^{1+} \\cup \\dotso \\cup S^{K+}$, which should just be the set of all sample indexes. Can you please make your notation more clear as to what you actually mean here?\n\nIn line 166 please clarify the notation $\\mathrm{mid}$.\n\nTypos/vocabulary:\n- line 83: it should be $\\mathbb R^{S \\times A}$\n- line 85: \"objective\"\n- line 124: \"belongings\" -> \"assignment\"\n- line 227: \"miscellaneous\"\n- line 228: \"symmetry-breaking\"", "limitations": "not applicable", "rating": 7, "confidence": 4, "soundness": 4, "presentation": 3, "contribution": 3, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "The paper makes a strong technical contribution to a computationally challenging problem. The writing is clear and concise, the numerical results are significant. Overall it's a good paper with only few flaws.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "7SQskK5FYm", "reviewer": "Reviewer_htGe", "summary": "The K-Medoids, as a classic clustering method, have many real applications. To fix the pool scalability of classic branch and bound (BB), this work proposed a method with reduced-space BB, in which the lower bound is provided by Lagrangian relaxation. Combined with several tailored bound tightening techniques, the method was tested in 28 datasets in the experiments section, including two datasets with one million samples.", "questions": "It would be better if the author can emphasize the significance of their upper/lower bound method.", "limitations": "No potential negative societal impact. ", "rating": 5, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "I do recognize the huge sets of experiments and the effort. Just saying, I have a feeling that: Compared to machine learning, this works probably more like the FOCS or STOC type of paper for theoretical computing.\n\nStrengths：\nThe idea is intuitive and makes sense.\n\nThere is an experiments-side contribution since the author compares the state-of-art baselines. The author provided experiments on 28 datasets, including two datasets with one million samples. The experiment on this data size is impressive.\n\nSeveral tailored bound tightening techniques needs effort as an engineering solution, especially when facing large-scale dataset.\n\nWeaknesses：\nTheoretical side, the focus should be on how we select the proper upper and lower bound in the BB. In this work, the upper bounding methods are fine but not surprising. And, as they mentioned, the tight lower bound from the work of Cornuejols. \n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "1AGf341KU2B", "reviewer": "Reviewer_FDa1", "summary": "This paper uses a branch-and-bound (BB) method to improve Cornuejols work on Lagrangian relaxation for the K-Medoids clustering problem. Through bounding the regions of the medoids, search space is significantly reduced, leading to notable computation speed improvements. A convergence proof is also provided. To obtain the lower and upper bound, a couple of tricks are used: 1. Other than the Lagrangian based lower bound, the author proposes a simple basic lower bound, to further speed up the process. 2. The author proposes heuristic methods to obtain the upper bound; 3. Different bound tightening methods (probing, FBBT) are applied before updating the search space. The proposed approach can run on millions of data samples within reasonable time. \n", "questions": "It would be better to introduce more on BB method in this setting, e.g. what are the BB tree and node? How are they constructed?", "limitations": "NA", "rating": 5, "confidence": 3, "soundness": 3, "presentation": 2, "contribution": 2, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\n1. The combination of brand-and-bounding with Lagrangian relaxation, on the K-Medoids problem is relatively novel. This method makes the Lagrangian based K-Medoids approach much more applicable in real-world datasets that have more than thousands of samples.\n2. The proposed approach is trivially parallelizable, which makes the method suitable for large scale problems.\n3. Theoretical convergence proof is provided.\n\nWeaknesses:\nThe biggest weakness of the paper is regarding the practical value of the proposed method. It might interest the optimization community due to its convergence guarantee and various bound-shrinking tricks, but as a clustering method, the proposed approach is significantly more complicated than the counterparts like K-Means:\n1. Though vanilla K-Means does not converge to global optimal, simple modifications like K-Means++ could alleviate the problem and often performs well in practice. The paper lacks a comparison (in terms of speed and quality) with such a simple baseline (K-Means or K-Means style methods).\n2. Other than the centroid-based approaches, there are many other kinds of clustering, e.g. density-based, hierarchical clustering, etc. They often don't need the ball-shape assumption and work better in real world datasets. The paper lacks comparison with these different approaches, especially the state-of-the-art ones like [Sarfraz, et al, Efficient Parameter-free Clustering Using First Neighbor Relations].\n3.  In experiments, the performance are evaluated using the achieved optimal gap and number of solved nodes. Standard clustering metrics include ARI and NMI, which provide a sense of clustering quality. The gap is useful in theoretical settings but not very meaningful in applications.\n\n------\n\nScore update from borderline reject to borderline accept, after reviewing the authors' responses. ", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}]

Author rebuttal:

Decision: Accept