Heuristic Methods for Large Centroid Clustering Problems.
This article presents new heuristic methods for solving a class of hard centroid clustering problems including the p-median, the sum-of-squares clustering and the multi-source Weber problems. Centroid clustering is to partition a set of entities into a given number of subsets and to find the location of a centre for each subset in such a way that a dissimilarity measure between the entities and the centres is minimized. The first method proposed is a candidate list search that produces good solutions in a short amount of time if the number of centres in the problem is not too large. The second method is a general local optimization approach that finds very good solutions. The third method is designed for problems with a large number of centres&semi; it decomposes the problem into subproblems that are solved independently. Numerical results show that these methods are efficientdozens of best solutions known to problem instances of the literature have been improvedand fast, handling problem instances with more than 85,000 entities and 15,000 centresmuch larger than those solved in the literature. The expected complexity of these new procedures is discussed and shown to be comparable to that of an existing method which is known to be very fast.
INTRODUCTION
.
analysis is to partition a set of entities into subsets, or clusters, such that the subsets are
homogeneous and separated one another, considering measurements describing the entities. This
problem already preoccupies Aristotle and appears in many practical applications. For instance, it
has been studied by naturalists in the XVIII th century for classifying living species. In this paper, we
propose new efficient methods for centroid clustering problems. More precisely, we are going to
apply our methods to problems of the following type: given n entities e i with weights w
it is searched p centres c j (j = 1, ., p) minimizing , where
measures the dissimilarity between e i and c j . However, the methods are very general and
may be applied to other problems or objective functions.
If the entities are described by their co-ordinates in IR m , d(e i , c j ) is typically the distance or the
square of the distance between e i and c j . In the last case, the problem is the well known
sum-of-squares clustering (SSC) (see e.g. Ward (1963), Edwards and Cavalli-Sforza (1965), Jancey
(1966), MacQueen(1967)). There are many commercial softwares that implement approximation
procedures for this hard problem. For instance, the popular S-Plus statistical analysis software incorporates
the k-means iterative relocation algorithm of Hartigan (1975) to try to improve the quality of
1. A former version of the article was entitled "Heuristic methods for large multi-source Weber problems".
given clusters. For exact algorithms for SSC, see e. g. Koontz, Narendra and Fukunaga (1975) and
Diehr (1985).
In case: 1) the space is IR 2 , i. e. the Euclidean plane, 2) the centres can be placed everywhere in
the dissimilarity measure is the Euclidean distance, the problem is called the multi-source
Weber problem (MWP). This problem occurs in many practical applications, such as the placement
of warehouses, emitter antennas, public facilities, airports, emergency services, etc. See e. g. Saaty
(1972), Dokmeci (1977), Fleischmann and Paraschis (1988), Bhaskaran (1992) and Lentnek,
MacPerson and Phillips (1993) that describe practical applications that need to solve MWPs with up
to more than 1700 entities and 160 centres. For exact methods solving the MWP, see e. g. Rosing
(1992) and Krau (1997). For a unified comparison of numerous approximation algorithms, see Brim-
berg et al. (1997).
In case dissimilarities between entities are given by an arbitrary n - n matrix and the centres can
be placed on the entities only, the problem is called the p-median problem (PMP). The last is a
well-known NP-hard problem, see e. g. Hakimi (1965), ReVelle and Swain (1970), Mirchandani and
Francis (1990) and Daskin (1995). For exact methods solving the PMP, see e. g. Erlenkotter (1978),
Rosing, ReVelle and Rosing-Vogelaar (1979), Beasley (1985) and Hanjoul and Peeters (1985). For
an introduction to location theory and clustering see also Gordon (1981), Sp-th (1985), Wesolowsky
(1993).
The new methods presented in this paper, candidate list search (CLS), local optimization (LOPT)
and decomposition/recombination (DEC), have been successfully applied to SSC, MWP and PMP,
but they can be extended to solve other problems. For example, the CLS and LOPT methods can be
applied to any location-allocation problems as soon as two appropriate procedures are available: the
first one for allocating entities to centres and the second one for optimally locating a centre, given the
entities allocated to it. For SSC, MWP or PMP, the allocation procedure simply consists in finding
the nearest centre to each entity. For other problems, this procedure must be more elaborated (e. g. if
there is a constraint limiting the sum of the weights of the entities allocated to a centre).
The LOPT method proceeds by local optimization of subproblems. This is a general optimization
method that can be applied to problems not directly related to clustering. For example,
multi-depot vehicle routing problems can be approached along its lines: the depots being identified
with the centres and the optimization of sub-problems being a procedure customized for solving relatively
small multi-depot vehicle routing problems. Also, the approach of Taillard (1993) for solving
large vehicle routing problems can be viewed as a special application of LOPT where the centres are
identified with the centres of gravity of the vehicle tours and the optimization procedure being an
efficient taboo search solving vehicle routing problems with few vehicles.
In order to remain relatively concise, we are going to present applications of our methods for
PMP, SSC and MWP only, but with a special attention to the under studied MWP. Indeed, while the
MWP by itself does not embrace all of the problem features found in some practical applications,
this model can be very useful, especially for real applications dealing with many thousands of enti-
3ties. In Figure 1, we show the decomposition into 23 clusters of a very irregular problem built on real
data, involving 2863 cities of Switzerland. The large black disks are the centres while the small disks
are the cities (or entities). Cities allocated to the same centres have the same colour. In this figure, we
have also added the federal frontiers and the lakes. Politically, Switzerland is composed of 23 states
physically, it is composed of extremely thickly populated regions (Plateau) and regions
without cities (Alps, lakes). We see in Figure 1 that the positions of the centres are sensible (no centres
are located outside Switzerland or on a mountain or in a lake) and that the decomposition generally
respects the natural barriers (spaces without cities); there are very few entities that are separated
from their centre by a chain of mountains 2 .

Figure

1 has to be compared with Figure 2 showing the decomposition of Switzerland into the
same number of clusters obtained by solving a PMP with dissimilarity measure being the true shortest
paths (the road network having more than 30000 connections). We see in this figure that the PMP
solution is very similar to the MWP one (21 centres are placed almost at the same position; the main
difference is that there are less entities allocated to a centre located on the other border of a lake).
However, solving this PMP is time consuming: the computation of the shortest paths matrix took
2. The expert can even identify a number of Swiss Cantons in this figure. There are however differences that could be
appropriate for solving political problems, such as the union of the South part of Jura to the Canton of Jura, the separation
of the German-speaking part of Valais or the union of the small primitive Cantons.

Figure

1: Decomposition of Switzerland into 23 clusters by solving a multi-source Weber Problem.
100 times longer than finding a very good MWP solution. Therefore, solving an MWP in a first
phase before attacking the true problem (as exemplified by a PMP or a multi-depot vehicle routing
problem) can be pertinent, even with an irregular, real problem.
Since the clustering problems treated in this paper are difficult, they can be solved exactly for
instances of moderate size only. For solving larger instances, as often arise in practice (see the 6800
entities, 380000 network nodes instance of Hikada and Okano (1997)), it is appropriate to use heuristic
methods. However, most of the methods of the literature present the same disadvantage of a
large increase of the computing time as the number of centres increases and, simultaneously, a
decrease in the quality of the solutions produced. The aim of this paper is to show that it is possible
to partition a problem with a large number of centres into subproblems that are much smaller, in
order to benefit from the advantages of the existing methods for small problems while rapidly producing
solutions of good quality to the original problem.
The article is structured as follows: in Section 2, we present in detail the alternate location-allocation
(ALT) procedure used as a subprocedure of our candidate list search (CLS), showing how it
can be implemented efficiently. ALT was first proposed by Cooper (1963) for the MWP. However, it
can be generalised for any location-allocation problem as soon as a location procedure and an allocation
procedure are available. In this section, we also present CLS, our basic procedure for solving the
subproblems generated by partition methods. In Section 3, we present two partition methods for
large problems. The first one, LOPT, can be viewed either as a generalization of the ALT procedure

Figure

2: Decomposition of Switzerland into 23 clusters obtained by solving a PMP with true shortest paths.
or as a restricted CLS for the post-optimization of a given solution. The second decomposition
method, DEC, splits a large problem into independent subproblems and the solutions of these sub-problems
are optimally mixed together to create a solution to the original problem. Section 4 analyses
the computational performances of the methods proposed.
2. BASIC PROCEDURES ALT AND CLS.
The procedures ALT and CLS are used as subprocedures in the decomposition methods we pro-
pose. Referring to the paper of Cooper (1963) is not sufficient to understand the procedure ALT well,
since certain details of this algorithm are not discussed in the original paper and the choices made for
implementing the procedure can have a profound impact on its effectiveness. Moreover, we have
adapted this procedure to accelerate its execution.
Generalized ALT procedure .
The iterative location-allocation procedure of Cooper (1963) may be sketched as follows:
Cooper has designed this algorithm for the MWP. In this case, the location procedure can be implemented
using a procedure like those of Weiszfeld (1937). For the SSC, the centre of gravity of the
entities is the optimum location of the centre. For the PMP, the optimum location of a centre can be
obtained by enumerating all possible location for the centre. The allocation procedure is very simple
for SSC, PMP and MWP: each entity is allocated to its nearest centre. For other problems, this procedure
can be more difficult to implement.
Two steps of this algorithm have to be discussed: the choice of the initial solution at step 1, and the
repositioning of centres that are not used at step 2a. For the choice of an initial solution, many variants
have been tested:
Position the centres on p randomly elected entities; the probability of choosing an entity being
proportional to its weight.
Choose the position of the centres one by one, by trying to position them on an entity and by
electing the position that minimizes the objective function.
The first variant takes into account the structure of the problem, i. e. the geographical and
weighting spread of the entities. It produces relatively good initial solutions, especially for problems
with non uniform weights.
Input: Set of entities with weight and dissimilarity measure,
problem specific allocation and location procedures.
Choose an initial position for each centre.
Repeat the following steps while the location of the centres
varies:
2a) Allocate the entities given the centre locations.
2b) Given the allocation made at step 2a, locate each centre
optimally.
Algorithm 1:Locate-allocate procedure of Cooper, 1963, (ALT).
The second variant induces the ALT procedure to produce the best solutions on the average but
its computing time is high: for each of the p centres, O(n) positions have to be tried, and for each of
these positions, one has to verify whether each entity is serviced by the new position. This implies a
procedure that operates in O(n 2 -p) time, while the other variant can be done much faster. To reduce
the complexity of this variant and to make it non deterministic 3 , we adopt the following O(n-p)
greedy procedure in the spirit of those of Dyer and Frieze (1985):
After having repositioned the centres at step 2b of the ALT procedure, it may happen that the allocation
of the next iteration, at step 2a, does not use all the centres. The unused centres can be relocated
to improve the current solution. We have adopted the following policy:
Determine the centre that contributes most to the objective function and place an unused centre
on its most distant entity; re-allocate the entities and repeat this as long as unused centres exist.
Starting with a very bad initial solution (O(p) centres that are not used), this re-location policy
could lead to a O(p 2 -n) procedure. However, our initial solution generator (as well as our CLS procedure
presented below) furnish solutions to the ALT procedure that contain an unused centre only
exceptionally (for the MWP, we have observed somewhat less than one occurrence in 1000, even for
a large number of centres). So, the re-location policy has almost no influence on the solution quality,
if one starts with a "good" initial solution as we do. Mladenovic and Brimberg (1996) have shown
that the re-location policy can have a substantial effect on MWP solution quality if one starts with
"bad" initial solutions.
Complexity of ALT for PMP, SSC and MWP.
First, let us introduce a new complexity notation: In the remaining of the paper, let -(.) denote
an empirically estimated complexity, while O(.) denotes the standard worst case complexity. For
example, both quick sort and bubble sort algorithms operate in O(n 2 ) time. In practice however, it is
3. In the context we use ALT, it is more interesting to have a non deterministic procedure. First, it may happen
that ALT is called many times for solving the same (sub-)problem. With a non deterministic procedure, it is
avoided to repeat exactely the same work. Then, let us mention that only non deterministic procedure can solve
NP-hard problems in polynomial time if P - NP. Therefore, our personal view is to consider non deterministic
procedure potentially more interesting than deterministic ones, even if there is no theory supporting this for the
moment.
Input: Set of entities with weight and dissimilarity measure.
Choose an entity at random and place a centre on this entity.
Allocate all entities to this centre and compute their weighted
dissimilarities.
3a) Find the entity that is the farthest from a centre (weighted
dissimilarities) and place the k th centre at that entity's
location.
3b) For is allocated to a centre farther
than centre k:
Allocate entity i to centre k and update its weighted dissimilarity

Algorithm 2:Initial solution generator.
observed that quick sort has an -(n-log(n)) behaviour while bubble sort has an -(n 2 ) behaviour
(Rapin, 1983) 4 . There are also algorithms for which the theoretical worst case complexity is not
established. However, observing the average computing times by executing an algorithm on many
instances can provide a good idea of its complexity in practice. The advantage of this notation is to
make a distinction between practice and theory. Indeed, it is common to read that the complexity of
quick sort is O(n-log(n)), which is not true, formally. Moreover, the "^" notation is often used by
statisticians for estimated values.
The complexity of the ALT procedure can be estimated as follows. The complexity of Step 2a
(allocation of the entities to a centre) is O(p-n). Indeed, for the problems under consideration one has
to allocate each entity to its nearest centre. For large values of p, this step can be substantially accelerated
by observing that only the centres that have moved from one iteration to the next can modify
the allocation previously made. (Compares the computing times of old and new ALT implementations
in Table 3.)
Step 2b can be performed in O(n) for the SSC. Indeed, each entity contributes only once in the
computation of the position of each centre (independently from the number of centres). For the
MWP, the optimum location can be found with a Weiszfeld-like procedure (1937) that repeats an
unknown number of gradient steps. We have arbitrarily limited this number to 30. So, in our imple-
mentation, Step 2b has a complexity of O(n). For small values of p, the computing time of this step
dominates. For the PMP, let us suppose that O(n/p) entities are allocated to each centre (this is reasonable
if the problem is relatively regular). 5 For each centre, one has to scan O(n/p) possible locations
and the evaluation of one position can be performed in O(n/p). So, the total complexity of
Step 2b is locating the p centres.
Since p is bounded by n, the global complexity of steps 2a and 2b is bounded by O(n 2 ) for SSC,
MWP and PMP. Now, we have to estimate the number of repetitions of Loop 2 which is unknown.
However, in practice, we have observed that the number of iterations seems to be polynomial in n
and p. Therefore, we will use an -(p a -n b ) estimation of the overall complexity of our implementation
of the ALT procedure. In this study, we are mostly interested in instances with large values of p,
so, we have considered instances with n/5 - p - n/3 for evaluating the a and b values for the various
clustering problems. For the SSC and MWP, we have considered about 7000 instances uniformly
generated with up to 9400 entities. For the PMP, we have considered about 38000 runs of the ALT
procedure. The PMP instances were based on the 40 different distance matrices proposed by Beasley
(1985). The number of entities for these instances ranges from 100 to 900.
4.More precisely, such a behaviour can be mathematically proven. In that case, we propose to follow the usual
notation in statistics and to write for an expected running time derived from a mathematical analysis. Therefore
it can be written that the complexity of quick sort is .
5.Without this assumption, the complexity is higher; with stronger assumptions (e. g. Euclidean distances), a
lower complexity can be derived.
O .
O
log
For the SSC, we have estimated a @ 0.83 and b @ 1.19; for the PMP the estimation is a @ 0.70
and b @ 1.23 and for the MWP a @ 0.85 and b @ 1.34. So, if p grows linearly with n, the estimated
complexity of the ALT procedure is not far from -(n 2 ) for all these problem types. The memory
requirement is O(n) for the SSC and MWP and O(n 2 ) for the PMP i. e. equivalent to the data size.
. 2 . Candidate list search (CLS) .
CLS is based on a greedy procedure that randomly perturbs a solution that is locally optimal
according to the ALT procedure. Then, ALT is applied to the perturbed solution and the resulting
solution is accepted only if it is better than the initial one, otherwise one returns to the initial solu-
tion. The perturbation of a solution consists in eliminating a centre and in adding another one,
located on an entity. The process can be repeated until all pairs entities/centres have been scanned.
This greedy procedure finds very good solutions: In Table 1, we report the quality of the solutions
found when applied to the 40 PMP instances of Beasley (1985). These instances have been solved
exactly and the quality of a solution is given in per cent above the optimum value. The greedy procedure
was executed 20 times for each instance. For 8 instances each run found the global optimum
and all instances but one were optimally solved at least once.
For MWP instances with 50 (respectively 287) entities we observed that the greedy procedure
finds a global optimum in more than 60% (respectively 40%) of the cases. For the SSC, we succeeded
in improving all the best solutions known to 16 instances with 1060 entities and 10 to 160
centres (See Table 4).
For the p-median problem, this type of perturbation has been used for a long time (c. f. Goodchild
and Noronha, 1983, Whitaker, 1983, Glover, 1990, Vo-, 1996, Rolland, Schilling and Current,
1997); in this case Glover proposes an efficient way to evaluate the cost of eliminating a centre: during
the allocation phase, the second closest centre is memorized - this can be done without increasing
the complexity. However, evaluating the decrease of the cost due to the opening of a centre on an
entity takes a time proportional to n. Therefore, finding the best possible perturbation has a complexity
of O(n 2 -p), without considering the application of the ALT procedure.
This complexity is too high for large instances thus we make use of a candidate list strategy
scheme proposed by Glover (1990) for implementing a probabilistic perturbation mechanism. The
idea is to identify the centre to close by a non deterministic but systematic approach. The entity associated
with an open centre is also randomly chosen, but its weighted distance from its previously

Table

1: Quality of the greedy procedure for Beasley's PMP instances (% above optimum).
allocated centre must be higher than the average. The process is repeated for a number q of itera-
tions, specified by the user. Algorithm 3 presents CLS into details.
The most time consuming part of this algorithm is step 3e, i. e. the application of the ALT procedure
to the perturbed solution. As seen above, we can estimate the complexity of this step as
-(p a -n b ). Therefore, the complexity of CLS is -(q-p a -n b ). From now on, we write CLS(q) the
improvement of a given solution with q iterations of the CLS procedure.
3. DECOMPOSITION METHODS.
In this section, we propose two decomposition methods for solving problems with a large
number p of centres. The complexity of these methods is not higher than the ALT procedure while
producing solutions of much higher quality. The first decomposition technique, LOPT, starts with
any solution with p centres and improves it by considering a series of subproblems involving r < p
centres and the entities allocated to them. The subproblems are solved by our CLS algorithm. This
method can be viewed as a local search defined on a very large neighbourhood involving up to r centres
re-locations at a time. Another point of view is to consider this procedure as a generalization of
ALT. Indeed, a solution produced by ALT is locally optimal if we consider any subset of entities
allocated to a single centre: every entity is serviced by the nearest centre, and the centres are optimally
positioned for the subset of entities they are servicing. Our procedure produces a solution that
is sub-optimal (since the subproblems are solved in a heuristic way and since we do not consider all
subsets of r centres) for subsets of entities allocated to r centres. A third point of view is to consider
LOPT as a CLS procedure with a much smaller list of candidate moves regarding to the CLS presented
above.
The second decomposition method, DEC, partitions the problem into t smaller subproblems.
These subproblems are then solved with our CLS for various numbers of centres. A solution to the
(location of the p centres), parameter
q.
Generate p, a random permutation of the elements {1, ., p} and
-, a random permutation of the elements
the distance of the most (weighted) distant
entity.
3c) While the weighted distance from entity - i to the nearest
centre is lower than (d
3d) Close centre p j and open a new one located at entity - i to
obtain a perturbed solution s
with ALT to obtain s k ''
new random permutation p.
Algorithm 3:Candidate list search (CLS).
initial problem is then found by combining solutions of the subproblems. To decompose the initial
problem, we solve an intermediate problem with t centres with our CLS procedure. Each set of entities
allocated to a centre of the intermediate problem is considered as an independent subproblem.
3 . 1 . Local optimization (LOPT) .
The basic idea of LOPT is to select a centre, a few of its closest centres and the set of entities
allocated to them to create a subproblem. We try to improve the solution of this subproblem with
CLS. If an improved solution is found, then all the selected centres are inserted in a candidate list C,
otherwise the first centre used for creating the subproblem is removed from C. Initially, all the centres
are in C and the process stops when C is empty. LOPT has two parameters: r, the number of centres
of the subproblems and s, the number of iterations of each call to CLS. Algorithm 4 presents
more formally the LOPT method.
Complexity of LOPT.
To estimate the complexity of LOPT we make two assumptions. First we assume that O(n/p)
entities are allocated to each centre (this hypothesis is reasonable if the problem instance is relatively
uniform) and second that loop 3 is repeated -(p g -n l ) times. Empirically, we have observed that g is
less than 1 and l is close to 0 (see Table 8); we estimate that the value of g is about 0.9 and l is about
for the LOPT parameters we have chosen and for the MWP). Then, the complexity of LOPT can
be established as follows:
Steps 3a and 3d have a complexity of O(p); step 3b has a complexity of O(r-p); step 3c solves a
problem with r centres and O(r-n/p) entities, this leads to a complexity of -(s-r a -(r-n/p) b ). This leads
to a total complexity of -(r-p are fixed and if p grows linearly
with n, the complexity of the LOPT procedure is therefore -(n l This complexity seems
to be similar to that of the ALT procedure. In practice, step 3c of the LOPT procedure takes most of
the computing time, even if steps 3b has a higher expected complexity for extremely large p. Indeed,
for fixed n, we have always observed that the computing time diminishes as p increases, even for p
larger than 10000 (see Tables 4 to 8). From now on, we denote by LOPT(r, s) the version of the
LOPT procedure using parameters r and s. The memory requirement of the LOPT procedure is O(n).
initial position of the p centres, parameters r and s.
While C -, repeat the following steps:
3a) Randomly select a centre i - C.
3b) Let R be the subset of the r closest centres to i (i - R).
3c) Consider the subproblem constructed with the entities allocated
to the centres of R and optimize this subproblem with
r centres with CLS(s).
3d) If no improved solution has been found at step 3c, set
else set
Algorithm 4: Local optimization procedure (LOPT).
3 . 2 . Decomposition algorithm (DEC) .
LOPT optimizes the position of a given number of centres dynamically, but it is also possible to
proceed to a static decomposition of the entities, and solve these subproblems with a variable
number of centres. A solution to the complete problem may be found by choosing the right number
of centres for each subproblem. Naturally, the total number of centres must be limited to p. This
re-composition may be performed efficiently and optimally with dynamic programming.
The crucial phase of the algorithm is the first decomposition: if the subproblems created do not
have the right structure, it is impossible to obtain a good solution at the end. The more irregular the
problem is (i. e. where the entities are not uniformly distributed, or their weights differ widely), the
more delicate its decomposition is. For partitioning the problem, we use our CLS procedure applied
to the same set of entities but with a number t < p of centres.
The subproblems created may have very different sizes: a subproblem may consist of just a few
entities with very high weights or it may comprise a large number of close entities. Thus it could be
difficult to evaluate the number of centres to be assigned to a subproblem. Let n i be the
number of entities of subproblem i. Suppose that subproblem i is solved with
and let f ij be the value of the objective function when solving subproblem i with j centres. To
build a solution to the initial problem, we have to find
This problem is a kind of knapsack and may be reformulated as:
Thus, the problem can be decomposed and solved recursively by dynamic programming in
O(t-p) time. This procedure can also produce all the solutions with t, t centres in O(t-n)
time. Such a feature can be very useful when we want to solve a problem for which the number of
centres is unknown and must be determined, as for example when there is an opening cost for each
centre (the opening cost has just to be added in the f ij values).
However, solving each subproblem with 1, ., n i centres is time consuming. If the problem is
relatively uniform, one can expect that the optimum number of centres found by dynamic programming
is not far from p/t for all subproblems. So, we propose to first solve the subproblems for only
three different numbers of centres: -p/t - 1-p/t- and -p/t + 1-. These are solved with one less
minimize
such that j i p
minimize
minimize
such that j i
(respectively one more) centre when the optimum number of centres determined by dynamic programming
is exactly the lower (respectively higher) number for which a solution was computed.
Algorithm 5 presents our DEC procedure in details.
Complexity of DEC .
For analysing the complexity of DEC, we make the following assumptions: First, each subproblem
has O(n/t) entities, second, each subproblem is assigned O(p/t) centres and third, the number of
repetitions of loop 5 is a constant (i. e. the total number of subproblems solved with CLS in steps 3
and 5a is in O(t) ). These assumptions are empirically verified if the problem instances are relatively
uniform (see Table 8). With these assumptions, the complexity of DEC can be established as
follows: Step 1 is in -(u-t a -n b ); steps 3 and 5a can be performed in -(v-t 1 - a - b -p a -n b ); finally, the
complexity of dynamic programming, in steps 4 and 5b is O(t-p). The overall complexity of DEC
strongly depends on the parameter t. As shown in the next section, the quality of solutions produced
by CLS slightly diminishes as the number of centres increases. We therefore seek to reduce the
number of centres in the auxiliary problem and in subproblems as much as possible. For this pur-
pose, we have chosen . The overall complexity of our implementation of DEC is
are constant and p grows linearly with n, the complexity
is lower than the ALT procedure. The
memory requirement is O(n 3/2 ). DEC requires more memory than CLS and LOPT, but the increase is
not too high and we have succeeded in implementing all the algorithms on a personal workstation.
From now on, we note DEC(u, v) the use of the DEC procedure with parameters u
and v.
Input: Set of entities with dissimilarity measure.
Solve an auxiliary problem with t centres with CLS(u).
2) The subsets of entities allocated to the same centre form t
independent subproblems.
For each subproblem do: Solve subproblem i with CLS(v) with
centres and update the
associated.
Find a collection j 1 *, ., j t * of optimum number of centres to
attribute to each sub-problem with dynamic programming.
While
5a) For all (i,
using CLS(v) and update the f iji associated.
5b) Find a new collection j 1 *, ., j t * of optimum number of centres
to attribute to each subproblem with dynamic programming

Algorithm 5:The decomposition procedure (DEC).
4. NUMERICAL RESULTS.
4 . 1 . Test problems .
For the numerical results presented in this section, we consider six sets composed of 654, 1060,
2863, 3038, 14051 and 85900 entities respectively. The 2863 entities set is built on real data: the
entities are the cities of Switzerland and the weight of each city is the number of inhabitants. This set
is denoted CH2863.
9 130936.12 160 178764596.3 170 260281.77 170 1890823.0 5000 254125361
14 84807.669 250 126652250.1 300 184832.94 300 1404028.4 10000 166535699
50 29338.011 700 41923352.49 800 94301.618 800 824127.49
70 21465.436 900 28112316.53 1000 78458.720 1000 725300.72
90 17514.423 2000 475580.22
100 16083.535 U1060 2500 409677.92

Table

2: Number of centres and best solution values of the MWP instances.
The other sets correspond to the travelling salesman problems that can be found under the names
of P654, U1060, Pcb3038 Brd14051 and Pla85900 in the TSPLIB compiled by Reinelt (1995). For
these sets, all entities are weighted to one and the dissimilarity between two entities is the Euclidean
distance (for PMP and MWP) or the square of the Euclidean distance (for the SSC). From these six
sets of entities, we have constructed a large collection of instances by varying p. In Table 2, we give
the values of p we have considered for each set, and the best MWP solution known associated with
each p.
All the best solutions known have been found during the elaboration of methods presented in
this paper; some have been reported earlier in Hansen, Mladenovic and Taillard (1996) or in Brim-
berg et al. (1997) for P654 and U1060. For P654, we were able to find the same best solution values
reported by Brimberg et al. for p - 60, and to find better values for p > 60; for U1060, we succeeded
in improving all the best solution values with the exception of we got the same value.
The best solutions published in Brimberg et al. were obtained by considering more than 20 different
methods, and running each of them 10 times. This last reference also reports the optimum solution
values of smaller problem instances with 50 and 287 entities. We were able to find all these optimum
solution values with our CLS method. So, we conjecture that many of the solution values given in

Table

2 for the smallest set of entities are optimal. For the larger sets, we think that small improvements
can be obtained.
The aim of Table 2 is to provide new MWP instances and to assert the absolute quality of our
methods: Indeed we think that providing the relative quality (measured in per cent over the best solution
value of Table 2) allows comparisons to be made more easily than providing absolute solution
values. Sometimes, the best solutions known have been found by using sets of parameters for which
results are not reported in this paper and it would be difficult to estimate the effort needed to obtain
each best solution known. Consequently, we do not provide computing times in this table.
Our algorithms are implemented in C++ and run on a Silicon Graphics (SG) 195MHz workstation
with R10000 processor. In order to make fair comparisons with algorithms implemented by
other authors and executed on a different machine, we have sometimes used another computer,
clearly indicated in the tables that follow. It was not possible to report exhaustive numerical results
due to the large number of problem instances (160), problem types (PMP, MWP or SSC) and methods
(CLS, DEC and LOPT). We try to report representative results in a condensed form. However,
let us mention that the conclusions we draw for a given method for a problem type are generally
valid for another problem type.
4 . 2 . ALT and CLS .
First, we want to show the efficiency of our CLS algorithm by comparing it to the results produced
by one of the best methods at the present time for the MWP: MWPM, an algorithm that first
solves exactly a p-median before re-locating optimally the centres in the continuous plane. This
method is due to Cooper (1963) but has been forgotten for a long time before Hansen, Mladenovic
and Taillard (1996) show that in fact, it is one of the most robust for small and medium size MWPs
(see also Brimberg et al., 1997). We do not consider methods such as those of Bongartz, Calamai and
Conn (1994) which are too slow and produce too poor solutions or those of Chen (1983) or Murtagh
and Niwattisyawong (1982) which are not competitive according to Bongartz et al. Also, we do not
compare our results with the HACA algorithm of Moreno, Rodr-gez and Jim-nez (1990) for two
reasons: First the complexity of HACA is O(p 2 n) and requires an O(p 2 ) memory, i. e. O(n 3 ) in time
and O(n 2 ) in memory if are clearly higher than those of our methods. Second,
produces solutions that are not as good as MWPM. Indeed, HACA first builds a heuristic
solution to the p-median instance associated to the MWP and then applies the ALT procedure to the
p-median solution. The reader is referred to Brimberg et al., 1997 for a unified comparison of a large
range of heuristic methods for the MWP.
To show the effects of the improvements of the ALT procedure proposed in this paper, we provide
the best solutions obtained over 100 repetitions of an old version of ALT that starts with different
initial solutions; this method is denoted MALT(100). The results for MALT(100) and MWPM
originate from Hansen, Mladenovic and Taillard (1996). In Table 3, we give the solution quality
(measured in per cent above the solution value given in Table 2) of MWPM, MALT(100), CLS(100)
and CLS(1000) and their respective computing times (seconds on Sun
P654. The computing time of CLS(1000) is roughly 10 times that of CLS(100). We have averaged
all these results for 10 independent runs of the algorithm. Where the 10 runs of CLS(100) find solutions
values identical to those given in Table 2, we provide in brackets the number of iterations
Quality (% above best known) Computing time [s. Sun Sparc 10]
7.3

Table

3: Comparisons of CLS(100) and CLS(1000) with MWPM and MALT(100) for MWP instances P654.
required by the worst run of CLS out of 10 to find the best solution known. From this table, we can
conclude that:
- The new ALT procedure runs 6 to 9 times faster than the old one (both MALT(100) and
call 100 times an ALT procedure, the old one for MALT, the new one for CLS).
provides much better solutions than MALT(100).
- CLS(1000) provides better solutions than MWPM, in a much shorter computing time.
As p grows, the solution quality of all the algorithms diminishes.
4 . 3 . Decomposition methods DEC and LOPT.
As LOPT requires an initial solution in input, we indicate the performances of LOPT when
applied to the solution produced by the DEC procedure.

Table

4 compares CLS(1000), DEC(20, 50), LOPT(10, 50) and 3 VNS variants (due to Hansen
and Mladenovic (1997)) for SSC instances built on entities set U1060. This table gives: The best
solution value known (found with our methods), the solution quality of the methods (per cent over
best known; VNS results originate from Hansen and Mladenovic), and their respective computing
times (seconds on Sun workstation). The computing time of VNS1 and VNS2 is 150 seconds
for all instances. VNS3 corresponds to the best over ten executions of VNS2; therefore, its
computing time is 1500 seconds. It is shown in Hansen and Mladenovic that all VNS variants are
more efficient than other methods of the literature, such as the k-means algorithm of Hartigan
(1975). For LOPT, we do not take into consideration the computing time of DEC to obtain the initial
solution. From this table, we can conclude:
- For small values of p, CLS provides better solutions than DEC, VNSs.
- For the largest values of p, DEC produces fairly good solutions and their quality seems not to
decrease as p increases.
value known
Quality (% above best known) Computing time [s. Sparc 10]
CLS DEC LOPT VNS1 VNS2 VNS3 CLS DEC LOPT
50 255509536.2 0.33 7.54 0.45 30.65 1.97 0.54 114.8 10.8 44.1
90 110456793.7 1.04 7.45 0.49 46.08 1.52 0.78 122.4 12.0 25.2
100 96330296.40 1.14 7.29 0.44 44.51 2.23 1.06 125.2 12.3 24.7

Table

4: Comparison of CLS(1000), DEC(20, 50), LOPT(10, 50) and various VNSs for SSC instances U1060.
- The solution quality of LOPT is always very good and seems to be somewhat correlated with the
initial solution quality (obtained here with DEC).
- Unexpectedly, the computing time of LOPT and DEC diminishes as p increases; this is undoubtedly
due to the small number of entities of U1060.
produces better solutions than VNSs in a much lower computation
time.

Table

5 show the effect of the parameters of DEC and LOPT by confronting DEC(20, 50),
DEC(20, 200), LOPT(10, 50) (starting with the solution obtained by DEC(20, 50)) and
(starting with the solution obtained by DEC(20, 200)). This table provides the solution
quality and the computing times (seconds on SG) for the MWP instance CH2863; all the results
are averaged over 10 runs. From this table, we can conclude:
- For small values of p, the quality of DEC(20, 200) is slightly better than DEC(20, 50) but the
computing times are much higher.
Starting with solutions of similar quality, LOPT(10, 50) and LOPT(10, 200) produces solutions
of similar quality but the computing time of LOPT(10, 200) is much higher.
- For larger values of p, the quality of DEC(20, 50) slightly decreases but the quality of
remains almost constant.
- The computing times of DEC and LOPT diminishes as p increases.
greatly improves the solution quality obtained by DEC.
- The methods seems to be very robust since they provide good results for instances with a very
irregular distribution of dissimilarities.
Quality [%] Computation time [s. on SG]
20, 50 20, 200 10, 50 10, 200 20, 50 20, 200 10, 50 10, 200
100 3.4 3.2 0.28 0.20
28 160 56 183
170 4.4 4.1 0.24 0.14 22 126 43 163
190 4.8 4.3 0.34 0.20 24 132 45 157
43 149
300 5.3 4.3 0.47 0.20
400 5.1 4.0 0.63 0.38 17 103 36 108
500 4.8 3.4 1.02 0.39
700 5.8 3.7 1.48 0.52 19 106
900 6.6 4.0 1.59 0.71
1000 7.1 4.7 2.18 1.10 20 108 33 69

Table

5: Quality and computing time of DEC and LOPT for different parameter settings for MWP instance CH2863
In tables 6 and 7, we compare DEC + LOPT to a fast variant of VNS, called RVNS, for SSC and
PMP instances built on entities set Pcb3038. RVNS results originate from Hansen and Mladenovic
(1997). We have adapted the LOPT parameters in order to get comparable computation times. For all
SSC, PMP and MWP instances, we succeeded in improving the best solutions published in this last
reference. In Table 6, we can see that DEC + LOPT is able to find better solutions than RVNS in
shorter computation times. For the PMP, RVNS seems to be faster than DEC and LOPT for the
smallest number of centres. However, let us mention that our implementation derives directly from
the MWP one and is not optimized for the PMP. For example we do not compute the distances only
once at the beginning of the execution and store them in a (very large) matrix. For large number of
centres, DEC + LOPT is again faster and better than RVNS.
In

Table

8, we provide computational results for our methods DEC(20, 50), DEC(20, 200) and
LOPT(7, 50) (applied to the solution obtained with DEC(20, 200)) for MWP instances Brd14051
and Pla85900. We give the following data in this table: The number n of entities, the number p of
centres, the solution quality obtained by DEC and LOPT (per cent over best known), the respective
computing times (seconds on SG), the proportion of sub-problems solved by DEC and LOPT. For
DEC, this proportion corresponds to the number of subproblems solved divided by . For
LOPT, this proportion corresponds to the number of subproblems solved divided by p. The results
are averaged for 5 runs for Brd14051 and the methods were executed only once for Pla85900. From

Table

8, we can conclude that:
- The solution quality provided by DEC slightly decreases as p increases; this is due mainly to the
decrease in the solution quality provided by the CLS procedure when solving the subproblems.
Quality [%] Time [s. Sparc 10]
200 21885997.1 2.44 5.55 0.90 159.9 48 44
300 13290304.8 2.50 7.22 1.44 229.3
400 9362179.2 3.35 7.56 1.70 165.0 43 27
500 7102678.4 2.85 7.47 1.73 204.4 43 25

Table

Comparison of DEC(20, 50), LOPT(6, 40) and RVNS for SSC instance Pcb3038.
Quality [%] Time [s. Sparc 10]
200 238432.02 1.23 4.12 0.74 107.6 106 187
500 135467.85 0.88 4.01 0.71 209.7 59 81

Table

7: Comparison of DEC(20, 50), LOPT(7, 50) and RVNS for PMP instance Pcb3038.
- The computing times of DEC and LOPT diminishes as p increases. However, we can observe an
increase in DEC computation times - as predicted by the complexity analysis - only for very
large values of p. For LOPT we cannot observe such an increase, meaning that solving the sub-problems
takes more time than finding close centres for generating the subproblems.
- The solution quality of LOPT is very good, generally well below 1% over the value of the best
solution known.
- The proportion of subproblems solved by LOPT diminishes as p increases, showing that g is
smaller than 1.
- The proportion of subproblems solved by DEC seems to be constant, as assumed in the complexity
analysis of section 3.2.
5. CONCLUSIONS
In this article we have proposed three new methods for heuristically and rapidly solving centroid
clustering problems. First, we propose CLS, a candidate list search that rapidly produces good solutions
to problems with a moderate number p of centres. Second, we propose LOPT, a procedure that
locally optimizes the quality of a given solution. This method notably reduces the gap between the
initial solution and the best solution known. The third method proposed, DEC, is based on decomposing
the initial problem into subproblems. DEC and LOPT are well adapted to solve very large
Quality [%] Time [s. SG] Proportion
(7, 50)100 2.4 2.38 0.39 458 1931 3109 4.4 4.4 4.4
200 1.9 2.03 0.30 336 1379 1632 4.3 4.6 3.3
300 2.2 2.08 0.30 252 1073 1055 5.0 5.0 2.6
400 2.2 2.04 0.26 214 915 885 4.9 5.2 2.4
500 2.5 2.12 0.27 195 909 838 4.9 5.8 2.5
2.4 1.99 0.26 184 799 707 5.1 5.8 2.3
700 2.5 1.93 0.23 177 829 632 5.7 6.8 2.2
800 2.5 1.97 0.24 149 760 555 4.8 6.8 2.1
900 2.6 2.00 0.26 146 736 505 6.2 7.3 2.1
1000 3.0 2.17 0.31 129 667 437 4.6 6.6 2.0
2000 4.0 2.81 0.59 88 379 227 4.7 4.8 1.8
3000 4.7 3.65 1.15 82 347 173 4.6 4.9 1.8
5000 4.4 3.59 1.28 73 309 123 4.3 4.5 1.51000 1.78 1.53 0.09 3557 9415 7634 4.2 5.2 2.9
1500 1.97 1.70 0.17 3149 7885 5343 4.3 5.7 2.8
2000 1.81 1.46 0.12 2819 6923 4750 4.1 5.1 2.5
3000 1.74 1.30 0.10 2405 5959 4532 4.1 5.2 2.6
5000 1.87 1.41 0.00 2597 5214 3423 4.1 4.3 2.0
7000 2.03 1.50 0.13 2276 4770 2526 4.3 4.3 1.9
8000 2.07 1.53 0.03 2681 4685 2344 4.1 4.6 2.0
9000 2.55 1.67 0.16 2796 4658 1992 4.1 4.3 1.8
10000 2.78 1.80 0.16 2629 4863 1813 5.1 4.5 1.8
15000 3.71 2.61 0.58 3144 5242 1552 5.3 6.0 1.8

Table

8: Computational results for DEC and LOPT for MWP instances Brd14051 and Pla85900.
problems since their computing time increases more slowly with the number of entities than that of
other methods in the literature. These methods can solve problems whose size is many order of magnitude
larger than the problems treated up to now. Despite its speed, they produce solutions of good
quality. The expected complexity of these procedures are given and experimentally verified on very
large problem instances.
In fact, LOPT is a general optimization method that can be considered as a new meta-heuristic.
Indeed it can be adapted for solving any large optimization problem that can be decomposed into
independent sub-problems. LOPT has been shown to be very efficient for centroid clustering problems
and vehicle routing problems. Future works should consider to apply LOPT to other combinatorial
optimization problems.
The success of the methods presented in this paper could be explained as follows: solving problems
with a very limited number of centres (e. g. below 15) is generally an easy task. Thanks to the
use of an adequate neighbourhood, the CLS method allows problems up to 50-70 centres to be
treated in a satisfactory way. DEC treats the problem at a high level and is able to determine the general
structure of good solutions involving a large number of centres. Starting with a solution that has
a good structure, LOPT is able to find very good solutions using a very simple improving approach.
Therefore, it is interesting to remark that a very simple improving scheme can lead to a very efficient
method if an initial solution with a good structure can be identified and an efficient neighbourhood is
used. Indeed, the quality of the solutions obtained by the DEC + LOPT method rival what one would
expect from a more elaborate meta-heuristic such as a genetic algorithm, taboo search or simulated
annealing. The use of inadequate neighbourhood structures can explain the poor performances of
previous implementation of such meta-heuristics. In summary, we can say that our methods open
new horizons in the solution of large and hard clustering problems.

Acknowledgements

The author would like to thank Fred Glover for suggesting numerous improvements, Ken Rosing
and Nicki Schraudolf for their constructive comments. This research was supported by the Swiss
National Science Foundation, project number 21-45653.95.
6.



--R

"A Note on Solving Large p-Median Problems"
"Identification of transshipment center locations"
"A projection method for l p norm location-allocation prob- lems"
"Improvements and Comparison of Heuristics for Solving the Multisource Weber Problem"
"Solution of minisum and minimax location-allocation problems with Euclidean distances"
"Location-allocation problems"
Network and Discrete Location
"Evaluation of a branch and bound algorithm for clustering"
"A quantitative model to plan regional health facility systems"
"A Simple Heuristic for the p-Centre Problem"
"A Method for Cluster Analisis"
"A dual-based procedure for uncapacitated facility location"
"Solving a large scale districting problem: a case report"
"Location-Allocation for Small Computers"
Classification: Methods for the Exploratory Analysis of Multivariate Data
"Tabu Search for the p-Median Problem"
"Optimum distribution of switching centers in a communication network and some related graph theoretic problems"
"An approach for solving a real-world facility location problem using digital map"
"A comparison of two dual-based procedures for solving the p-median prob- lem"
"Heuristic solution of the multisource Weber problem as a p-median problem"
"An introduction to Variable Neighborhood Search"
Clustering Algorithms
"Multidimensional Group Analysis"
"A Branch and bound clustering algorithm"
Extensions du probl-me de Weber
"Optimum producer services location"
"Some Methods for Classification and Analysis of Multivariate Observations"
Discrete Location Theory

"Heuristic cluster algorithm for multiple facility location-allocation problem"
"An efficient method for the multi-depot location-allocation problem"
Cours d'informatique g-n-rale
"TSPLIB95"
"An efficient tabu search procedure for the p-Median Prob- lem"
"An Optimal Method for Solving the (Generalized) Multi-Weber Porblem"
"The p-Median and its Linear Programming Relaxation: An Approach to Large Problems"
"Optimum positions for airports"
Dissection and Analysis (Theory
"Parallel iterative search methods for vehicle routing problems"
"A reverse elimination approach for the p-median problem"
"Hierarchical Grouping to Optimize an Objective Function"
"Sur le point pour lequel la somme des distances de n points donn-s est minimum"
"The Weber problem: history and perspectives"
"A fast algorithm for the greedy interchange for large-scale clustering and median location problems"
--TR

--CTR
Hongzhong Jia , Fernando Ordez , Maged M. Dessouky, Solution approaches for facility location of medical supplies for large-scale emergencies, Computers and Industrial Engineering, v.52 n.2, p.257-276, March, 2007
Mauricio G. C. Resende , Renato F. Werneck, A Hybrid Heuristic for the p-Median Problem, Journal of Heuristics, v.10 n.1, p.59-88, January 2004
Teodor Gabriel Crainic , Michel Gendreau , Pierre Hansen , Nenad Mladenovi, Cooperative Parallel Variable Neighborhood Search for the p-Median, Journal of Heuristics, v.10 n.3, p.293-314, May 2004
multi-source Weber problem;sum-of-squares clustering;p-median;clustering;location-allocation