# Using $k$-medoids for distributed approximate similarity search with
arbitrary distances

Elena Garcia-Morato
CETINIA-DSLAB
Universidad Rey Juan Carlos
Madrid, Spain
<EMAIL_ADDRESS>
& Maria Jesús Algar
CETINIA-DSLAB
Universidad Rey Juan Carlos
Madrid, Spain
<EMAIL_ADDRESS>
& Cesar Alfaro
CETINIA-DSLAB
Universidad Rey Juan Carlos
Madrid, Spain
<EMAIL_ADDRESS>
& Felipe Ortega
CETINIA-DSLAB
Universidad Rey Juan Carlos
Madrid, Spain
<EMAIL_ADDRESS>
& Javier Gomez
CETINIA-DSLAB
Universidad Rey Juan Carlos
Madrid, Spain
<EMAIL_ADDRESS>
& Javier M. Moguerza
CETINIA-DSLAB
Universidad Rey Juan Carlos
Madrid, Spain
<EMAIL_ADDRESS>
Corresponding author. This article is currently under peer-review. All
datasets are either publicly available or can be completely reproduced
following the guidance provided in the associated references.

###### Abstract

This paper presents GMASK, a general algorithm for distributed approximate
similarity search that accepts any arbitrary distance function. GMASK requires
a clustering algorithm that induces Voronoi regions in a dataset and returns a
representative element for each region. Then, it creates a multilevel indexing
structure suitable for large datasets with high dimensionality and sparsity,
usually stored in distributed systems. Many similarity search algorithms rely
on $k$-means, typically associated with the Euclidean distance, which is
inappropriate for specific problems. Instead, in this work we implement GMASK
using $k$-medoids to make it compatible with any distance and a wider range of
problems. Experimental results verify the applicability of this method with
real datasets, improving the performance of alternative algorithms for
approximate similarity search. In addition, results confirm existing
intuitions regarding the advantages of using certain instances of the
Minkowski distance in high-dimensional datasets.

_Keywords_ Approximate similarity search $\cdot$ k-medoids $\cdot$ arbitrary
distance function

##  1 Introduction

The goal of similarity search is to find elements from a dataset, represented
by feature vectors in a multidimensional space, similar to a given query
object according to a particular similarity or dissimilarity measure. Numerous
data indexing methods, collectively known as access methods, have been
proposed to speed up this task. Nevertheless, the results remain strongly
conditioned by the problems inherent to the representation of datasets in
high-dimensionality and high-sparsity spaces. These algorithms are frequently
conceived to work with data having spatial components or represented in a
vector space with few dimensions. However, they become unusable in high-
dimensional spaces, where many features characterise each element due to the
Curse of Dimensionality [1]. As a result, the distance between elements tends
to become increasingly uniform, making it challenging to differentiate between
them based solely on their distances or when the number of dimensions
significantly exceeds the quantity of data points. Moreover, several studies
indicate that in datasets with over 20 dimensions, exact similarity search
methods rarely outperform exhaustive scanning.

Owing to this constraint, novel approaches have emerged in recent years to
swiftly identify enough related items within a reasonable time frame compared
to linear scanning, often referred to as approximate similarity search. One of
the newest approaches is MASK (Multilevel Approximate Similarity search with
k-means) [2], a novel indexing method that involves an unconventional
application of the $k$-means partitioning algorithm [3, 4] to create a
multilevel index structure for approximate similarity search.

Traditionally, the input to $k$-means is always a fixed number of $k$
prototype points or centroids. Then, following an iterative process, the
dataset is divided into $k$ groups (Voronoi cells), assigning each point to
its closest centroid according to a given distance function, so that each
centroid represents all of the elements of that cluster. In contrast,
leveraging current powerful computing capabilities, $k$-means can be forced to
generate a high number of centroids that inundate the dataset. Then, these
centroids become a new dataset to be grouped at the upper level, repeating
this process until the complete construction of the search tree is achieved.
This bottom-up creation of the multilevel index enables the distribution of
this task among different computing nodes, which also store distinct data
partitions.

Although it employs an innovative approach that capitalises on the benefits
offered by contemporary computing systems and distributed platforms, the
performance of this indexing method remains significantly impacted by the
inherent constraints of the $k$-means algorithm. The choice of the distance
function significantly impacts the outcome of $k$-means clustering, as it
determines how the similarity between data points is calculated and can lead
to different outcomes in terms of the shape and separation of the clusters.

Until now, experiments conducted with $k$-means have relied on the Euclidean
distance to calculate relative distances between points, following the
conventional inclination to associate the concept of distance with that used
in geospatial contexts. Nonetheless, relying on the Euclidean distance as the
default option in indexing structures and algorithms when working with high-
dimensional data, despite its widespread use, can lead to sub-optimal results.
This is because, in high-dimensional spaces, the discriminating power of the
$L_{p}$ norm decreases for higher values of $p$ due to their specific features
(such as dispersion, outliers, and concentration of points around certain
values), leading to poor differentiation between proximate and distant
elements.

This paper introduces GMASK (General Multilevel Approximate Similarity Search
with $k$-medoids), a novel algorithm that facilitates the incorporation of
arbitrary distances by employing data partitioning algorithms that induce
Voronoi regions in a dataset and yield a representative element, exemplified
by $k$-medoids. The simple introduction of $k$-medoids turns GMASK into a
versatile algorithm capable of adapting to the characteristics of any dataset
(high dispersion, high dimensionality, outliers, etc.) and seek for the
optimal results. This approach allows using multiple $\ell_{p}$ norms or even
other distance functions, making it compatible with a broader spectrum of
problems. To comprehensively evaluate its performance, an extensive series of
experiments will be conducted across datasets varying in type, size,
dimensionality, and sparsity.

This paper adheres to the following structure. Section “Background and related
work” reviews previous research related to similarity search. Section
“Methodology” describes GMASK, emphasising the potential for enhancing the
algorithm’s performance by exploring the application of other distance
functions beyond the Euclidean one, like other Minkowski distances. The
findings of four empirical experiments carried out over real-world datasets
with quite diverse attributes are presented in Section “Experimental results”.
These tests compare GMASK performance with other well-known algorithms, such
as FLANN and PyNNdescent. Finally, Section “Conclusions and future work”
recaps the main conclusions of this work and suggests promising lines for
future work.

###  1.1 Contributions

The potential for achieving improved results, coupled with the lack of prior
empirical studies that validate this proposal and the need for learning
paradigms applicable to a wider range of dissimilarity measures [5], justifies
our efforts to explore alternative Minkowski distances for data indexing. This
endeavour is aimed at challenging the misconception we have identified and
fostering a deeper understanding of the significant impact of using
alternative distance metrics when working with high-dimensional data.

  * •

GMASK allows the use of different clustering algorithms for the bottom-top
construction of the index tree, such as $k$-medoids, a broader form of the
$k$-means algorithm that offers parameterisation with several distance
metrics, in particular Minkowski distances.

  * •

Evaluating the performance of different $L_{p}$ norms or Minkowski distances
for approximate $k$-nearest neighbour search. This assessment involves testing
GMASK against other state-of-the-art methods to illustrate that the selection
of a suitable distance function, while considering the features of the target
dataset, has a significant impact on performance when searching for reasonably
similar regions.

##  2 Background and related work

This section provides minimal background on the concepts and strategies
underlying similarity search methods, including the description of current
similarity search algorithms.

###  2.1 Similarity search

Similarity search, often referred to as proximity search [6, 7, 8], is a
cornerstone for applications in many different fields, such as databases [9],
information retrieval [10], computer vision [11] and bioinformatics [12, 13],
among others. It aims to find elements that are similar (close) to a given
query object, according to some measure of similarity or, conversely,
dissimilarity.

Thus, given a query object $q\in\mathbb{X}$, the similarity search tackles
some of the following scenarios, represented in Figure 1 [7]:

  * •

Point query: Finding elements within $\mathbb{X}$ that have exactly the same
feature values as $q$.

  * •

Range query: Identifying a subset of elements $\\{o_{i}\\}\in\mathbb{X}$ whose
feature values do not exceed a given similarity threshold $r$ around $q$ such
that $\forall\\{o_{i}\\},\delta(o_{i},q)\leq r$.

  * •

Nearest-neighbour query: Finding elements in $\mathbb{X}$ whose features are
the closest to the feature vector representing $q$. In this case, the goal
might be to retrieve the single element (set of $k$ elements) closest to $q$,
known as the nearest neighbour ($k$-nearest neighbours) and denoted by
$NN_{q}$ ($k-NN_{q}$). Identifying the $k$-nearest neighbours of a given
object is one of the most widely used classification methods in machine
learning.

Figure 1: Conceptual representation of several types of query searches.

Many techniques for data indexing, collectively known as access methods, can
speed up these queries. These methods create an index structure that
partitions the domain $\mathbb{X}$ into different regions based on a distance
function. Then, we can leverage these regions to solve similarity queries.
There is a wide range of access methods, including exact approaches,
approximate solutions and variants tailored for distributed computing systems.

###  2.2 Algorithms for similarity search

Exact similarity search methods identify the completely accurate result set
for a query by using the absolute positions of objects to group elements and
search for similar instances. These methods are categorised as Point Access
Methods (PAM) for point-based elements and Spatial Access Methods (SAM) for
elements with spatial dimensions, such as lines, polygons or higher-
dimensional polyhedra. The general acronym SAM is commonly used in the
research literature to refer to both categories of indexing strategies in
vector spaces [10]. These approaches are most effective with elements
represented in a relatively small number of dimensions, making Curse of
Dimensionality one of the biggest challenges in data science today.

Metric Access Methods (MAM) provide a general framework for similarity search
problems within metric spaces [6, 8], exploiting the triangle inequality to
partition the metric space into subspaces, effectively filtering out portions
of the dataset that cannot contain valid results. On the one hand, MAM exact
search can be expensive in terms of computational time and updating the index
structure due to dynamic data. On the other hand, MAM indexes have serious
problems in efficiently reducing the search for candidate elements in a query
when dealing with high-dimensional data. There are many MAM for queries in
metric spaces, some of which use clustering algorithms, such as $k$-means [4]
or $k$-medoids [14], to recursively partition the entire set of elements
creating a hierarchical structure (tree) of nested clusters.

Approximate indexing methods are solutions that strike a balance between
accuracy and speed, making them suitable for both vector and metric spaces
[15]. Although these methods may introduce some errors into the results, they
also provide a faster response time by calculating fewer distances. As a
result, these methods prove feasible in scenarios where exact methods cannot
provide a fast response, such as in high-dimensional problems or large
datasets, where exact methods default to a full scan. Some of the most widely
used exact algorithms are compared and summarised in [16], which also
introduces a comprehensive taxonomy to classify them according to relevant
characteristics.

Rather than computing the exact result set for a query, approximate similarity
search algorithms adopt various strategies (e.g., modifying the representation
space, reducing the number of pair distances to compute) to identify a reduced
set of candidate items that are likely to be the most similar. FastMap,
introduced in [17], is one of the first algorithms to be developed for
approximate similarity search in metric spaces. It is based on a method to
obtain a fast projection of elements in the original dataset into a user-
defined $k$-dimensional space, preserving the distances between pairs of
elements in order to achieve an approximate response through dimensionality
reduction. In this category of algorithms, MASK [2] emerged as a novel
indexing approach characterised by a multilevel design tailored for
approximate similarity search. It features an innovative utilisation of the
$k$-means partitioning algorithm to quickly reduce the number of checked
regions when searching for results.

###  2.3 Using different distance functions

The dissimilarity value captures the notion of closeness between two objects,
which is interpreted as a distance in a suitable space. A distance function
that avoids introducing any correlation between features, such as the
Minkowski distance [18] or $\ell_{p}$ norm, can quantify this value [17]. The
Manhattan or city block distance ($\ell_{1}$ norm), the Euclidean distance
($\ell_{2}$ norm), and the Chebyshev distance ($\ell_{\infty}$ norm) are
typical instances of this family. As many multidimensional indexing structures
were initially proposed in the context of spatial systems, many indexing
structures and algorithms have commonly relied on the Euclidean distance,
derived from its conventional use in two- or three-dimensional spatial
applications.

However, several negative consequences and challenges arise when dealing with
high dimensional problems, an effect widely known as the Curse of
Dimensionality. As the number of dimensions increases, the volume of the space
grows exponentially, and the density of points in the space decreases. This
sparsity poses a challenge in accurately measuring distances and similarities
between data points, since proximity becomes less meaningful. Consequently,
points can be very far apart, even though they are close according to a metric
such as Euclidean distance. This effect increases the computational complexity
of algorithms and operations performed in high-dimensional spaces, often
resulting in slower and less efficient computations. For this reason, several
authors have pointed out that, due to the Curse of Dimensionality, it is very
costly to find the exact nearest neighbour in high-dimensional Euclidean
spaces [19, 20].

As a result, although Euclidean spaces reflect our natural conception of
space, this natural assumption is challenging to extrapolate when dealing with
high-dimensional data. Moreover, traditional clustering algorithms, such as
$k$-means, often specify the Euclidean distance as the only appropriate metric
without considering other possibilities. This choice brings in a self-imposed
limitation that may produce results far from the best [21, 22, 5].

In recent years, several surveys [23, 24, 25] have addressed the evaluation of
various metrics for similarity searches across spaces of diverse topology, a
burgeoning topic offering insights into the most suitable similarity measures
for different problems and contributing to the ongoing development of new
methods for similarity search.

Therefore, we aim to introduce GMASK, an algorithm that, being based on a
generalisation of earlier concepts previously introduced by MASK, signifies a
substantial departure from the latter. GMASK enables the integration of
arbitrary distances through the use of data partitioning algorithms that
generate Voronoi regions and assign a prototype for each of these regions.
Specifically, the utilisation of the $k$-medoids algorithm allows for the
incorporation of multiple $\ell_{p}$ norms, facilitating parametrisation with
different distance metrics and the execution of validation experiments. Our
findings demonstrate the feasibility of employing diverse distances beyond the
Euclidean metric, highlighting the potential for achieving superior results in
data indexing, contingent upon the specific characteristics of the dataset.

##  3 Methodology

This section describes the methodology developed for data indexation using
arbitrary distance functions. First, it illustrates the GMASK algorithm as a
general approximate indexing method accepting any distance function.
Hereafter, it delves into the characteristics of Minkowski distances,
emphasising the inherent limitations in $k$-means that justify the
implementation of GMASK. Consequently, we justify adopting $k$-medoids, an
alternative clustering method that leverages the opportunities provided by
Minkowski distances to overcome previous restrictions derived from using the
Euclidean distance. The design of experimental tests is described in
Subsections “Performance evaluation”, “Benchmark algorithms for similarity
search”, and “Datasets”, including metrics to assess performance, state-of-
the-art algorithms for comparison, and the datasets selected to gauge
effectiveness.

###  3.1 GMASK

As mentioned above, there is a close connection between machine learning and
similarity search indexing. Firstly, the foundations for designing data access
methods rely on familiar concepts in machine learning classification, such as
dissimilarity spaces and distance functions. Secondly, clustering methods like
$k$-means or $k$-medoids are used in many of these algorithms to create
indexing structures for similarity search and improve search performance.
Besides, they are also affected by the same limitations, such as the Curse of
Dimensionality and the challenge of identifying the actual intrinsic
dimensionality governing some problems to make them solvable.

The typical choice in partitioning clustering methods is to find a low value
of $k$ that generates different groups in the dataset. GMASK follows an
unconventional approach for clustering, generating a high number of $k$
representative points for each level. The only requirement in GMASK is using a
clustering method that generates Voronoi regions and returns a representative
point for each of them. In this implementation, GMASK uses the $k$-medoids
clustering method, inducing a Voronoi tessellation and returning an element
representing each region. Thus, generating a large number of representative
points, GMASK summarises the underlying data points with finer detail. In the
same way, this indexing method applies to any dissimilarity space, including
the particular case of metric spaces or vector spaces. Thanks to the
$k$-medoids adaptation in GMASK, this algorithm presents an important
improvement over MASK: it is able to work with different distance functions.
This modification brings GMASK the capability of adaptation to the specific
characteristics of any dataset. In the “Experimental results” section, we show
that this unusual application of clustering has attractive properties and
yields superior performance, even with high-dimensional datasets residing in
feature spaces that are typically sparse.

At this point, it is essential to note that the goal of GMASK is not to
cluster data but to use the partitioning clustering method to construct a
multilevel index to find elements. With more powerful computational
infrastructures and large memory and storage capacity becoming prevalent, the
advantage of this method lies in that previous restrictions on the number of
clusters at each level become less relevant. In this new scenario, we can
exploit the effect of introducing a high number of prototype points into the
dataset, leveraging the capabilities of the computing infrastructure.

Figure 2 illustrates the design principles of GMASK. At the lowest level, we
have the data points in $\mathcal{X}$ to be indexed. The blue boxes in the
diagram represent different data partitions or groups that could be stored on
various computing nodes. The novelty lies in creating the multilevel structure
of $k$ medoids, following a bottom-up approach. The first level of medoids
(identified as $n_{1}$ in Figure 2) summarises the data points at the lowest
level so that points assigned to each cluster are represented by its prototype
(medoid). The second level (denoted $n_{2}$) corresponds to medoids
representing the prototypes computed at the first level, $n_{1}$. The same
strategy is repeated, creating additional prototype layers that summarise the
immediately lower group of medoids. The procedure stops when the number of
medoids at the top level of this index structure becomes manageable.

Figure 2: Conceptual representation of the updated multilevel structure
approach created by GMASK.

The process of constructing this multilevel index structure depends on two
initial parameters:

  * •

The group length quantifies the number of elements in each data partition at
the lowest level, configured according to the available computing power (and
the total number of nodes in distributed computing systems). For the baseline,
we assume that all partitions are of the same size, although this condition is
not strictly necessary.

  * •

The number of medoids calculated for each group. At this point, the question
of selecting the most optimal number of medoids has not yet been addressed,
but it is not very important. Contrary to the standard procedure, our approach
suggests that, for information indexing, one should use as many medoids as the
computational infrastructure can afford.

The clustering at the first level reduces the initial number of data points by
a particular proportion, the data summarisation ratio, which is determined by
the ratio between the group length and the number of medoids. This ratio also
determines the number of layers of the multilevel index and the stop condition
when the number of medoids at the top layer becomes lower than the group
length.

Each medoid is a prototypical point representing all the elements assigned to
its cluster. The number of points assigned to each cluster may vary depending
on their distance from the corresponding medoid. For this reason, regions with
higher data density will also tend to receive more medoids. By generating a
high number of medoids, we can capture finer details of the underlying data,
allowing us to index not only high-density regions but also outliers and
clusters of elements distinct from the main concentrations of points.
Algorithm 1 describes the construction of this multilevel index structure.

Algorithm 1 Multilevel index construction in GMASK

1:function Clustering($distFunc,nPrototypes,points$)

2: $\triangleright$ A clustering method returning Voronoi regions

3: $\triangleright$ and a prototype point representing each region

4:end function

5:

6:procedure GMASK($data,lengthGroup,nPrototypes$)

7: $ngroups\leftarrow lengthData/lengthGroup$;

8: $vector\leftarrow$ split($data,nGroups$); $\triangleright$ Divide dataset
into nGroups

9: $\triangleright$ Initialise lists $layerPoints$, $layerLabels$, and
variable $idLayer$

10: while $ngroups\geq 1$ do

11: $\triangleright$ Initialise groupsPoints, groupsLabels and points

12: for $idGroup=1,nGroups$ do

13: $points\leftarrow vector[idGroup]$;

14: $groupsPoints.add($clustering($distFunc,nPrototypes,points$));

15: $groupsLabels.add(clusters.labels)$;

16: end for

17: $layerPoints.add(groupsPoints)$;

18: $layerLabels.add(groupsLabels)$;

19: $ngroups\leftarrow length(groupsPoints)/lengthGroup$;

20: if $ngroups\geq 1$ then

21: $vector\leftarrow$ split($layerPoints,nGroups$);

22: end if

23: $idLayer\mathrel{+}=1$;

24: end while

25: $layers\leftarrow idLayer-1$;

26:end procedure

Once built according to the rules described above, the index is ready to
accept point queries, which systematically search for all points in the
original dataset, seeking a perfect match. To solve any query using this
multilevel index, we start at the top level of the medoids hierarchy, steering
the search according to the distance between the medoids at each level and the
target query object $q\in\mathcal{X}$.

###  3.2 Dissimilarity measures

A dissimilarity space [26] defines the rules for measuring the distance
between elements represented in a feature space. As discussed in Section
“Similarity search”, similarity captures the notion of closeness between two
objects. Euclidean vector spaces are widely applied in previous research, as
they reflect and intuitive interpretation of a space: vectors are elements
that can be added, multiplied by a scalar or projected onto each other.

Being $\mathbb{X}$ a set of elements, the function

$\delta:\mathbb{X}\times\mathbb{X}\to\mathbb{R}$ (1)

returns a real number $\delta(a,b)$ that shall be used to quantify the
difference between two or more elements contained on this set, denoting a
measure of distance.

Distance functions have the following properties:

  * $(p1)\quad\forall a,b\in\mathbb{X},\delta(a,b)\geq 0$ positiveness,

  * $(p2)\quad\forall a,b\in\mathbb{X},\delta(a,b)=\delta(b,a)$ symmetry,

  * $(p3)\quad\forall a\in\mathbb{X},\delta(a,a)=0$ reflexivity,

  * and in most cases

  * $(p4)\quad\forall a,b\in\mathbb{X},a\neq b\to\delta(a,b)>0$ strict positiveness.

If it also satisfies

  * $(p5)\quad\forall a,b,c\in\mathbb{X},\delta(a,b)\leq\delta(a,c)+\delta(c,b)$, commonly known as triangle inequality,

$\delta$ is indeed a metric, and the pair $(\mathbb{X},\delta)$ is called a
metric space [6]. If the distance does not satisfy the strict positiveness
property $(p4)$ then the space is called a pseudo-metric space, as well as if
it does not hold the symmetry property $(p2)$, we may have a quasi-metric
space. Nevertheless, there are some primitive spaces where these metric
constraints are not necessarily satisfied. We also use the word
’dissimilarity’ when referring to this difference or distance between elements
[18].

Normed and inner product spaces are special cases of metric spaces, where a
metric is defined either by a norm or an inner product. An Euclidean space is
a specific case of a Hilbert space. In turn, a Hilbert space is an inner
product space and a particular case of a Banach space, an example of normed
spaces in which a norm metric $\ell norm$ has been defined [5].

###  3.3 Minkowski distance

Minkowski distances or $\ell_{p}norms$ are distance measures commonly used to
compute the dissimilarity of objects described by numerical attributes.
However, dissimilarity can also be quantified using non-numerical attributes
through other alternatives, such as the edit distance widely used for text
attributes.

Let $\ell_{\infty}^{p},p\geq 1$ be a vector space of real sequences such that
$\sum_{i=1}^{\infty}|x_{i}|^{p}$ with the norm given by
$||x||_{p}(\sum_{i=1}^{\infty}|x_{i}|^{p})^{\frac{1}{p}}$. This norm induces
the $\ell norm$ or Minkowski metric $d_{p}$ and consequently $\ell_{p}^{m}$
and $\ell_{p}^{m}=(\mathbb{R}^{m},d_{p})$ are Banach spaces.

Thus, assuming $x,y\in\mathbb{R}^{m}$, the $\ell norm$ or Minkowski distance
is defined as follows:

$d_{p}(x,y)=(\sum_{i=1}^{m}|x_{i}-y_{i}|^{p})^{\frac{1}{p}},p\geq 1$ (2)

where depending on the $p$ value we can induce:

  * $d_{1}(x,y)=\sum_{i=1}^{m}|x_{i}-y_{i}|$, the Manhattan or city-block distance.

  * $d_{2}(x,y)=d_{E}(x,y)=(\sum_{i=1}^{m}|x_{i}-y_{i}|^{2})^{\frac{1}{2}}$, the Euclidean distance.

  * d${}_{\infty}(x,y)=d_{m}ax(x,y)=max_{i}|x_{i}-y_{i}|$, the max-norm or Chebyshev distance.

Minkowski distances with $p<1$ cannot be considered as proper metrics since
they do not hold the triangle inequality. Nonetheless, several authors have
shown that fractional distances can significantly improve the effectiveness of
standard clustering algorithms such as $k$-means [22].

###  3.4 Limitations of k-means and the Euclidean distance

In unsupervised machine learning with $k$-means, the user must determine the
appropriate number of clusters for a given problem, and a centroid represents
each cluster. Carefully choosing the number of centroids is crucial because
too many centroids would not effectively summarise the underlying data, and
too few centroids would result in grouping unrelated data points. However, the
locations provided to these initial centroids can greatly impact the final
outcome, making it a sensitive aspect.

To address this limitation, algorithms such as $k$-means++ [27] propose an
approach where the initial random centroids are spread more evenly, leading to
improved results. Nevertheless, although several more powerful $k$-means
implementations have tried to reduce the performance impact resulting from
sub-optimal initial centroid placement, only a few have explored the potential
for achieving better solutions by choosing a more appropriate distance.

Due to the nature of the $k$-means algorithm, the choice of distance function
can also exert a significant impact on $k$-means outcomes, as it determines
how similarity between data points is calculated and produces clusters with
different shapes and separations between each other. As it has been raised in
Section “Background and related work”, despite the widespread use of Euclidean
distance in indexing structures and algorithms, relying on it as the default
option when using the $k$-means algorithm can lead to poor results.

Previous research [22] suggests that inherent characteristics of high-
dimensional spaces, such as high dispersion, outliers or concentration of
points around specific values, hinder the discriminatory power of the
$\ell_{p}norm$ for higher values of $p$. This situation would provide undue
importance to the distant (sparse/noisy) dimensions and poorly discriminate
between further and nearer neighbours. Therefore, for a given high-dimensional
problem, it may be preferable to use Minkowski distances with lower values of
$p$, so the $\ell_{1}$ distance metric (Manhattan distance) would be the best
choice for high-dimensional applications, followed by the Euclidean distance
($\ell_{2}$), then the $\ell_{3}$ metric, up to the Chebyshev distance
($\ell_{\infty}$).

Previous algorithms for approximate similarity search, such as MASK [2], have
relied on the Euclidean distance to calculate relative dissimilarities between
points. This strategy follows the conventional inclination to associate
distance between points with the Euclidean metric, which aligns with the
distance concept commonly applied in geospatial problems. However, the
potential for achieving improved results, coupled with the lack of prior
empirical studies that validate this proposal and the need for learning
paradigms applicable to a broader range of dissimilarity measures [5], justify
our efforts to explore alternative Minkowski distances for data indexing. This
endeavour aims to foster a deeper understanding of the impact of using
alternative distance metrics with high-dimensional data.

Nevertheless, while it may seem theoretically plausible to employ $k$-means or
its enhanced version, $k$-means++, with several distance metrics, the
practical implementation reveals that the algorithm’s internal design is
intertwined with the Euclidean distance. This alignment relies on the
algorithm’s iterative optimisation process when it searches for centroids that
minimise the total squared distances of points within their respective
clusters[3]. The sum of the Squared Error (SSE), calculated as the sum of
squared Euclidean distances between each data point and the centroid of its
assigned cluster, underscores the algorithm’s inherent reliance on Euclidean
geometry for clustering optimisation [28].

As we have already demonstrated, assuming that widely adopted mechanisms are
the most suitable or the only viable option can lead to self-imposed
limitations that generate sub-optimal outcomes. Instead, our objective is to
steer clear of the misconception that certain implementations are the best or
the only conceivable choice simply because they are the most common choices.

At this juncture, GMASK can integrate an alternative clustering algorithm that
leverages the opportunities provided by diverse distance metrics for data
indexing. The $k$-means algorithm defines clusters based on centroids
minimising the sum of squared errors in distances from all other points in the
same cluster. In contrast, alternative methods like $k$-medoids calculate
clusters around the most representative data element, which is determined by
minimising the sum of pairwise dissimilarities among data points the same
cluster.

###  3.5 k-medoids as a generalisation of k-means

Several generalisations of the $k$-means algorithm have been proposed to
broaden the scope of applicable data types [29], reflecting the diversity
found in nature prototype points, such as subsets of objects, hyperplanes,
etc.

The first paper proposes a general criterion and its generalised $k$-means
method is [30]. This approach starts from a $n$ x $n$ dissimilarity matrix
$(d_{\ell}t)$ for $n$ factories $O=1,...,n$ in an industrial network where
$d_{\ell}t$ are the minimum road transportation costs between $\ell$ and $t$.
The goal is to partition the set of factories into $k$ classes
$(C_{1},...,C_{k})$ and to find a selection $Z=(z_{1},...,z_{k})$ of $k$
factories as ’supply points’ such that, when supplying all factories of the
class $C_{i}$ from the supply point $z_{i}\in O$, the overall transport costs
are minimised, where $d(\ell,z_{i})=d_{\ell},z_{i}$ means the dissimilarity
between the factory (object) $\ell$ and the factory (supply point) $z_{i}\in
O$. So the family $P$ of admissible ’prototypes’ consists of all singletons
from $O$ and the class prototypes, meaning ’most cheapest supply point’ in
$C_{i}$, are determined by:

$Q(C_{i},z):=\sum_{\ell\in C_{i}}d(\ell,z)\to\min_{z\in C_{i}}$ (3)

Many authors, including, among others, Diday [31, 32, 33], Bock [34] or
Govaert [35], have followed the generalised clustering approach in various
settings and numerous variations, obtaining a plethora of generalised
$k$-means algorithms.

Kaufmann and Rousseeuw [36] first propose the popular Partitioning Around
Medoids (PAM) algorithm, also referred to as $k$-medoids. In $k$-medoids, data
is modelled by using $k$ representative objects $m_{i}$ called medoids, that
serve as prototypes for the clusters (formerly called classes) which are
obtained by using the absolute error criterion (’total deviation’, TD). Since
the TD is defined as the sum of dissimilarities of each point $x_{j}\in C_{i}$
to the medoid $m_{i}$ of its cluster, the medoid of a cluster $C_{i}$ is in
fact the object with the smallest sum of dissimilarities (or equivalently, the
smallest average) to all other objects in the cluster.

$medoid(C_{i}):=argmin_{x_{i}\in C}\sum_{{x_{j}\in C}}d(x_{i},x_{j})$ (4)

This method is suitable for non-Euclidean distances and a wide range of data
since it does not force the dissimilarity measure to adhere strictly to metric
properties, and only requires a proximity measure between any pair of objects.
In addition, while $k$-means defines a prototype in terms of a centroid, which
is usually the mean of a group of points that may not correspond to an element
in the dataset, $k$-medoids defines a prototype in terms of a medoid, an
actual data point identified as the most representative one for a group of
points.

Regarding $k$-medoid implementations, Program PAM [36] consists of two
algorithms, BUILD to choose an initial clustering, and SWAP to improve the
clustering towards a local optimum, since finding the global optimum of the
$k$-medoids problem is, unfortunately, NP-hard. The algorithms require a
dissimilarity matrix, which requires $O(n^{2})$ memory and at least
$O(n^{2}d)$ time to be computed. On top of this, the process of updating
medoids involves evaluating the cost of swapping non-medoid points to ensure
that the sum of pairwise dissimilarities within clusters (rather than the sum
of squared Euclidean distances) is minimised, making this method more robust
to noise and outliers compared to the alternatives based on $k$-means.
Recently, FasterPAM [37] introduced a speedup for larger $k$, achieving an
$O(k)$-fold acceleration in the second SWAP phase of the algorithm, but still
providing the same results as the original PAM algorithm.

Originally FasterPAM was implemented in Java and published within the open-
source library ELKI [38]. However, to integrate it with GMASK implementation,
developed entirely in Python, we will employ the fast $k$-medoids Python
implementation [39] version provided by the $kmedoids$ package
111https://pypi.org/project/kmedoids. FasterPAM let us upgrade our method by
making it compatible with any distance function and a wider range of problems.

Experimental results verify the applicability of GMASK with real datasets,
improving the performance of alternative algorithms and providing better
results for approximate similarity search queries.

###  3.6 Performance evaluation

The main advantage of GMASK is providing a fast response to queries at the
expense of returning approximate results. Thus, a straightforward strategy to
assess the performance of this algorithm is to undertake an exhaustive search
of all elements in a given dataset $X$. Then, the proportion of data points
that were correctly retrieved can be determined. This is the common approach
to all evaluation experiments developed in this study.

A comprehensive series of experiments are conducted across datasets of diverse
types, sizes, dimensionality, and sparsity to evaluate the performance of
GMASK, which introduces $k$-medoids for the bottom-top construction of the
index tree. This novelty enables its parametrisation with several $\ell_{p}$
norms as outlined in “k-medoids as a generalisation of k-means”. These results
will be compared with those obtained using other state-of-the-art algorithms.
As mentioned above, our objective is to understand better the impact of
utilising distances other than the Euclidean metric. Through these
experiments, we aim to demonstrate that data indexing using non-Euclidean
distances is not only feasible but also has the potential to yield improved
outcomes.

###  3.7 Benchmark algorithms for similarity search

To evaluate our thesis, we compare GMASK performance in approximate nearest
neighbour search with respect to alternative MAM.

####  3.7.1 KD-Tree

The KD-Tree (K-Dimensional Tree) [40] is an exact similarity search method
used to find the nearest neighbours to a given query point. It constructs a
binary tree structure by recursively dividing the data based on the maximum
spread along different dimensions and the search algorithm efficiently narrows
down the search space by traversing the tree. In low-dimensional spaces, KD-
trees can achieve sublinear search times. However, as the dimensionality
increases, the performance of the KD-tree can degrade due to the Curse of
Dimensionality phenomenon.

As it performs an exact search, the solutions provided by this method are used
as a baseline to evaluate the performance of other approximate algorithms. We
shall compare the k-nearest neighbours obtained by each method against those
obtained using the KD-Tree to assess their accuracy.

####  3.7.2 FLANN

FLANN (Fast Library for Approximate Nearest Neighbours) [11] is an efficient
method widely employed for approximate similarity search in large-scale and
high-dimensional datasets, extensively used in machine learning, computer
vision, and data mining. FLANN employs techniques such as random projection
and hierarchical subdivision to construct an index structure that enables
faster search operations. It offers several index types, including KD-tree,
$k$-means tree, and hierarchical clustering tree, which are suitable for
diverse data types and query requirements. While FLANN has implementations in
several programming languages, we use the Python bindings provided by the
pyflann 222https://github.com/primetang/pyflann library.

####  3.7.3 PyNNDescent

The PyNNDescent algorithm is an approximate similarity search method [41]
implemented in Python 333https://pynndescent.readthedocs.io/en/latest/. It is
based on the NNDescent algorithm and used for efficient nearest neighbour
search in large-scale datasets. PyNNDescent constructs a graph-based index by
connecting each data point to its approximate nearest neighbours. It
iteratively improves the graph structure through neighbour descent steps,
employing techniques such as random initialisation, multi-tree searching, and
pruning. The algorithm also provides parameters to control the trade-off
between search accuracy and computational cost. PyNNDescent is particularly
useful for high-dimensional datasets and is widely used in machine learning,
data mining, and information retrieval tasks.

###  3.8 Datasets

Several experiments have been conducted to evaluate the performance of GMASK,
using four datasets diverse in nature and public access. Table 1 briefly
summarises their main features.

####  3.8.1 Municipalities

The Municipalities dataset is a two-dimensional geospatial dataset that
contains the coordinates (latitude and longitude) of each of the 8,130
municipalities that make up the Spanish territory, as depicted in Figure 3. It
has been extracted from one of the tables provided in the ’2020 Resumen
Descriptivo del Nomenclátor de Municipios y Entidades de Población’
(Descriptive Summary of the Nomenclature of Municipalities and Population
Entities), called MUNICIPIOS444https://ign.es/web/rcc-nomenclator-nacional.
This report is compiled annually by the Instituto Geográfico Nacional
(National Geographic Institute), a public institution under the Spanish
Ministry of Transport and Sustainable Mobility. It includes updated geographic
and administrative data on all Spanish municipalities, such as their name,
province, geographic location coordinates, perimeter or population, among
others.

Figure 3: Visualisation for dataset Municipalities

Although this is a two-dimensional dataset with a low dispersion, it is an
exceptionally intriguing case due to its intuitive representation and the
presence of outliers, such as the Canary Islands and the Balearic Islands. The
location of these islands deviates from the rest of the points significantly
and poses a challenge to any similarity search method applied to this dataset.

####  3.8.2 MNIST

The MNIST_784 dataset 555http://yann.lecun.com/exdb/mnist/ is a widely used
machine learning benchmark dataset whose simplicity makes it an ideal starting
point for algorithm development and evaluation. Firstly, it consists of 60,000
training images and 10,000 test images, a relatively small size that makes
possible to carry out computationally efficient experimentation and
prototyping. Secondly, it comprises grayscale images of handwritten digits
ranging from 0 to 9, making it visually straightforward and easy to work with.
Moreover, each image is 28x28 pixels in size, resulting in a total of 784
dimensions when flattened, allowing researchers to explore algorithms in the
context of high-dimensional data. Many well-established techniques and
algorithms have been evaluated on MNIST [42, 43], making it a reference point
for evaluating the performance of approximate similarity search methods [16,
44].

####  3.8.3 GLOVE

The third experiment involves the GLOVE_100
dataset666https://nlp.stanford.edu/data/glove.6B.zip, which contains one
million pre-trained Global Vectors for Word Representation (GloVe) embeddings
that capture semantic relationships between words. These 100-dimensional
embedding vectors were pre-trained on the combined Wikipedia 2014 + Gigaword
5th Edition corpora (6B tokens, 400K vocabulary) and cover a wide range of
English words, enabling comprehensive assessments of approximate similarity
search methods across diverse vocabulary and word relationships. Given their
widespread adoption as a benchmark in the natural language processing
community, the use of GloVe embeddings as a standardised dataset provides
valuable insights into the method effectiveness in retrieving similar items
based on their semantic similarity [16].

####  3.8.4 NYtimes

The NYtimes dataset, part of the Bags of Words dataset from the UCI
repository777https://archive.ics.uci.edu/dataset/164/bag+of+words, comprises a
collection of New York Times news articles represented as a bag of words. Each
document in the dataset is associated with a set of word occurrences, where
the dimensions represent unique words extracted from the articles. The dataset
is organised as a document–word matrix, where each row corresponds to a
document and each column corresponds to a word. The values in the matrix
indicate the frequency of each word occurring in the respective document.
Preprocessing steps include tokenization, removal of stopwords, and vocabulary
truncation, with only words occurring more than ten times retained.

The reduced version of the NYtimes employed on this work was generated by
building tf-idf descriptors from the bag-of-words version, and embedding them
into a lower dimensional space using the Johnson-Lindenstrauss Transform [45].

Dataset |  Label |  n |  Dimensionality | High Sparsity |  Data Type
---|---|---|---|---|---
Municipalities |  $\mathcal{X}_{muni}$ |  8,130 |  2 | No |  Geospatial
MNIST |  $\mathcal{X}_{MNIST}$ |  69,000 |  784 | Yes |  Image
GLOVE |  $\mathcal{X}_{GLOVE}$ |  1,000,000 |  100 | No |  Text
NYtimes |  $\mathcal{X}_{NYtimes}$ |  290,000 |  256 | No |  Text
Table 1: Descriptive traits characterising the datasets included in
experiments.

##  4 Experimental results

In order to show the versatility of GMASK, we present the analysis of four
different datasets: municipalities, MNIST, GLOVE and NYtimes. This study
consists of a neighbours search applying Euclidean distance, Manhattan
distance and Chebyshev distance. The main objective is to measure the accuracy
of the GMASK algorithm with different distances and to observe if any distance
is more adequate for some dataset, depending on its characteristics. For each
dataset, we compare the results obtained by GMASK with two other algorithms:
PYNN and FLANN. Both PYNN and GMASK are able to work with these three
distances, however, FLANN does not work with Chebyshev distance. Since MASK is
a specific case of GMASK only able to work with Euclidean distance, its
inclusion in these experiments has not been considered.

Our approach to evaluating GMASK performance is conceived to closely emulate
real-world conditions. The datasets used in this work are not artificially
generated, and no preprocessing has been applied. The transformation and
cleaning of the collected data is limited to dividing the original dataset
into two distinct subgroups, the training and testing datasets. We also check
that identical elements do not appear in both subgroups, thereby intensifying
the complexity in the search for $k$ approximate neighbours.

All experiments have been performed on a server equipped with 2 AMD EPYC 7451
microprocessors (24 cores/48 threads, 2.30 GHz, 64 MB L3 cache), 128 GB of
DDR4 RAM and an Intel D3-S4510 SSD (capacity 480 GB) for secondary storage.
The experimental results are presented in the following sections.

###  4.1 Municipalities

The municipalities dataset is an interesting case due to the geographical
distribution of points. In Figure 3, we can identify three main groups of
points: the largest region corresponds to the peninsular territory of Spain,
whereas the group of points located at the right of the map corresponds to the
Balearic Islands and the group of points located at the bottom left corner of
the map corresponds to the Canary Islands. The last two groups of points
increase the difficulty of the analysis because they are far from the rest of
the points, so we can consider them as outliers. Despite this dataset
containing additional information, we only consider two variables: latitude
and longitude.

Figure 4 shows the recall comparison between FLANN, PYNN and GMASK. We use a
heatmap to represent the recall value obtained for each algorithm with each
distance in the search of 5, 10 and 15 neighbours. The recall obtained in this
case for each method is very similar, and all of them depict good recall
values, over 84%. It is noticeable that FLANN recall decreases by about 10%
with the number of neighbours. Regardless of the distance used and the number
of neighbours intended to search, GMASK is the only method that provides a
100% recall in all tests.

Figure 4: Heatmap comparing the Recall (percentage of true nearest neighbours
found) obtained by each algorithm and distance in approximate nearest
neighbour search, for different number of neighbours, with dataset
$\mathcal{X}_{muni}$.

###  4.2 MNIST

Given its high dimensionality and sparsity, MNIST has become a widespread
benchmark in similarity search research. The following experiment evaluates
the recall obtained in searching neighbours in this dataset, applying
Manhattan, Euclidean and Chebyshev distance with FLANN, PYNN and GMASK
methods. Figure 5 presents the results of this analysis.

From this comparison, we can deduce that the results are highly comparable in
the recall when applying PYNN and GMASK methods with Euclidean and Manhattan
distances, respectively. Therefore, the combination of both distances and
these methods suits this case study. However, this is not the case for
Chebyshev distance since the recall is lower than 70%. There is also a
meaningful reduction in FLANN recall when the number of neighbours increases.
Regardless of the distance, GMASK yields a higher average recall than other
methods.

Figure 5: Heatmap comparing the Recall (percentage of true nearest neighbours
found) obtained by each algorithm and distance in approximate nearest
neighbour search, for different number of neighbours, with dataset
$\mathcal{X}_{MNIST}$.

###  4.3 GLOVE

Hereafter, the third experiment involves the analysis of the GLOVE dataset, a
high dimensional datasets that contains the largest number of elements among
the datasets included in these experiments.

Figure 6 depicts the recall results obtained for each method and distance.
Once again, the luck of efficiency is evident in FLANN´s results. The best
recall obtained with FLANN is 12% when searching for 5 neighbours applying
Euclidean distance. PYNN method provides good accuracy when Manhattan distance
is used. If we focus on PYNN method and compare the recall for the three
distances, we can see that Manhattan distance provides the best results and
Chebyshev distance provides the worst recalls. However, in the case of GMASK,
the recall is 100% per each distance and with independence of the number of
neighbours searched.

Figure 6: Heatmap comparing the Recall (percentage of true nearest neighbours
found) obtained by each algorithm and distance in approximate nearest
neighbour search, for different number of neighbours, with dataset
$\mathcal{X}_{GLOVE}$.

###  4.4 NYtimes

Finally, the last experiment replicates the same analysis using a dataset
characterised by high dimensionality, featuring a data-to-dimensionality ratio
lower than that of GLOVE but higher than MNIST: the NYtimes dataset.

Figure 7 illustrates the recall accomplished by FLANN, PYNN and GMASK. There
is a noticeable difference between GMASK and the other methods, since the
recall values for FLANN and PYNN are below 50%. However, GMASK consistently
outperforms the other methods regardless of the distance chosen. Once again,
FLANN exhibits the poorest performance.

Therefore, while the most remarkable insight obtained from this experiment is
the superior performance of the Chebyshev distance metric compared to others,
it is crucial to highlight that, in this case, GMASK nearly doubles the recall
provided by PYNN, achieving an 82% recall when searching for the 5 nearest
neighbours.

Through this experiment, we empirically demonstrate what was previously stated
theoretically in earlier sections of this manuscript: adopting the Euclidean
distance as the default option in indexing structures and algorithms is a
self-imposed restriction that can lead to sub-optimal results.

Figure 7: Heatmap comparing the Recall (percentage of true nearest neighbours
found) obtained by each algorithm and distance in approximate nearest
neighbour search, for different number of neighbours, with dataset
$\mathcal{X}_{NYtimes}$.

##  5 Discussion

Given the notable prevalence of Euclidean distance in similarity searches in
the previous literature, as well as the limitations detected in the widely
adopted $k$-means algorithm, which is strongly linked to this metric, there is
a growing need to explore the use of alternative metrics when addressing these
types of problems. Stemming from this need, GMASK makes the accommodation of
arbitrary distances feasible by integrating the $k$-medoids algorithm as a
clustering algorithm, allowing us to leverage the advantages derived from
using several distance metrics on high-dimensional datasets.

It is crucial to emphasise that the algorithm’s generalisation, achieved
through the integration of $k$-medoids, significantly broadens the range of
applicable distances while preserving the structural integrity and properties
of the original algorithm. Thus, GMASK maintains crucial characteristics,
including the desirable scalability property facilitated by the bottom-top
construction of the search tree, retaining a competitive edge over other
state-of-the-art algorithms that require centralising all data on a single
server.

As elaborated in Section “Experimental results”, our findings demonstrate the
efficiency of this new approach, thereby supporting the hypothesis that
alternative distances beyond the Euclidean metric may yield superior results
in solving nearest-neighbour queries, contingent upon the specific
characteristics of the dataset. Therefore, it becomes evident that the
outcomes achieved using the proposed method surpass those obtained through two
other widely adopted state-of-the-art algorithms, which are distinct from each
other and significantly differ from GMASK.

These findings would also be in line with previous publications [22]
suggesting that, for a given high-dimensional problem, using an $\ell_{p}$
norm with lower values of $p$ may be advantageous in assessing dissimilarity
between query points. Thus, while Aggarwal carries out a theoretical approach
using synthetic data, our research delves into this hypothesis by generalising
an algorithm that has been tested over a range of real-world datasets of
diverse nature. However, further empirical research is required to validate
this observed behaviour more precisely.

GMASK provides an avenue for future research to integrate distance functions
beyond Minkowski, thus providing many opportunities to extend this research
line. For instance, it would be interesting to explore the suitability of
distance functions for specific purposes, such as cosine similarity for
textual data or the Jaccard distance for categorical data, to address
particular problems.

Furthermore, although GMASK has yielded notably satisfactory experimental
results, there still exists potential for enhancing this implementation to
achieve further optimisation in computational methods, thereby increasing
overall performance. To this end, the proposal suggests exploring the
incorporation of a just-in-time (JIT) compiler for Python, such as Numba
888https://numba.pydata.org—employed by other algorithms such as PYNNDescent—,
or considering the adaptation and port of GMASK to a language renowned for its
strong performance in building scalable and concurrent systems like Go
999https://go.dev.

##  6 Conclusions and future work

Approximate similarity search is essential for many application areas where
large and heterogeneous datasets arise, often exhibiting high dimensionality
and sparsity. In many cases, these datasets must be stored in distributed
systems to become manageable, imposing additional restrictions that rule out
algorithms conceived for cases where data must reside on a single machine.

This paper presents the GMASK algorithm for approximate similarity search in
distributed systems that can accept arbitrary distance functions. One of the
main novelties is adopting the $k$-medoids (PAM) algorithm, which is not
bounded to the Euclidean distance, to achieve efficient solutions.
Furthermore, using FastPAM, a highly efficient implementation of the PAM
algorithm to solve the $k$-medoids problem, GMASK can improve its efficiency
with respect to other previous algorithms.

Experiments conducted on real-world datasets with different properties confirm
that GMASK, based on $k$-medoids, performs better in solving nearest-neighbour
queries than other frequently used algorithms, using other Minkowski
dissimilarity functions different from the Euclidean distance. These results
illustrate that indexing data with non-Euclidean distances is not only viable
but also holds the potential to enhance search outcomes substantially.
Overlooking alternative distances would represent a voluntary limitation that
may hinder us from achieving the optimal results.

GMASK lays out new interesting routes for further research. Once the
requirement of using the Euclidean distance has been lifted, additional
studies should be conducted to evaluate GMASK performance in other settings
where alternative distances can be applied, such as pseudo-metrics (Minkowski
distances with $p<1$) or non-metric dissimilarities. Besides, using arbitrary
distances will enable us to consider datasets with qualitative features
exclusively or a combination of quantitative and qualitative features. Many
state-of-the-art algorithms and standard dissimilarity measures cannot be
directly applied in these cases since they are incompatible with categorical
features.

## References

  * Böhm et al. [2001] Christian Böhm, Stefan Berchtold, and Daniel A. Keim.  Searching in High-Dimensional Spaces: Index Structures for Improving the Performance of Multimedia Databases.  _ACM Comput. Surv._ , 33(3):322–373, September 2001.  ISSN 0360-0300.  doi:10.1145/502807.502809.
  * Ortega et al. [2023] Felipe Ortega, Maria Jesus Algar, Isaac Martín de Diego, and Javier M. Moguerza.  Unconventional application of k-means for distributed approximate similarity search.  _Information Sciences_ , 619:208–234, 1 2023.  ISSN 0020-0255.  doi:10.1016/j.ins.2022.11.024.
  * MacQueen [1967] James MacQueen.  Some methods for classification and analysis of multivariate observations.  In _Proceedings of the Fifth Berkeley Symposium on Mathematical Statistics and Probability_ , volume 1: Statistics, pages 281–297, Berkeley, CA, USA, 1967. University of California Press.  URL https://projecteuclid.org/euclid.bsmsp/1200512992.
  * Lloyd [1982] Stuart P. Lloyd.  Least squares quantization in PCM.  _IEEE Trans. Inf. Theory_ , 28(2):129–136, 1982.  doi:10.1109/TIT.1982.1056489.
  * Pekalska and Duin [2005] Elzbieta Pekalska and Robert P. W. Duin.  _The dissimilarity representation for pattern recognition: foundations and applications_.  Number v. 64 in Series in machine perception and artificial intelligence. World Scientific, Hackensack, N.J, 2005.  ISBN 978-981-256-530-3.
  * Chávez et al. [2001] Edgar Chávez, Gonzalo Navarro, Ricardo Baeza-Yates, and José Luis Marroquín.  Searching in Metric Spaces.  _ACM Comput. Surv._ , 33(3):273–321, September 2001.  ISSN 0360-0300.  doi:10.1145/502807.502808.
  * Samet [2006] Hanan Samet.  _Foundations of Multidimensional and Metric Data Structures_.  The Morgan Kaufmann Series in Data Management Systems. Academic Press, San Francisco, CA, 2006.  ISBN 978-0-12-369446-1.
  * Zezula et al. [2006] Pavel Zezula, Giuseppe Amato, Vlastislav Dohnal, and Michal Batko.  _Similarity Search - The Metric Space Approach_ , volume 32 of _Advances in Database Systems_.  Springer, New York, NY, 2006.  ISBN 978-0-387-29146-8.  doi:10.1007/0-387-29151-2.
  * Agrawal et al. [1993] Rakesh Agrawal, Christos Faloutsos, and Arun N. Swami.  Efficient similarity search in sequence databases.  In David B. Lomet, editor, _Foundations of Data Organization and Algorithms, 4th International Conference, FODO’93, Chicago, Illinois, USA, October 13-15, 1993, Proceedings_ , volume 730 of _Lecture Notes in Computer Science_ , pages 69–84, Berlin, Heidelberg, 1993. Springer.  doi:10.1007/3-540-57301-1_5.
  * Baeza-Yates and Ribeiro-Neto [2011] Ricardo Baeza-Yates and Berthier Ribeiro-Neto.  _Modern Information Retrieval: The Concepts and Technology behind Search_.  Addison-Wesley Publishing Company, England, 2 edition, 2011.  ISBN 9780321416919.
  * Muja and Lowe [2014] Marius Muja and David G. Lowe.  Scalable Nearest Neighbor Algorithms for High Dimensional Data.  _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , 36(11):2227–2240, 2014.  doi:10.1109/TPAMI.2014.2321376.
  * Gish and States [1993] Warren Gish and David J States.  Identification of protein coding regions by database similarity search.  _Nature genetics_ , 3(3):266–272, 1993.
  * Tuncbag et al. [2011] Nurcan Tuncbag, Attila Gursoy, Ruth Nussinov, and Ozlem Keskin.  Predicting protein-protein interactions on a proteome scale by matching evolutionary and structural similarities at interfaces using PRISM.  _Nature Protocols_ , 6(9):1341–1354, 2011.  ISSN 1750-2799.  doi:10.1038/nprot.2011.367.
  * Kaufman and Rousseeuw [1990] Leonard Kaufman and Peter J. Rousseeuw.  _Finding Groups in Data: An Introduction to Cluster Analysis_.  Wiley Series in Probability and Statistics. John Wiley & Sons, Hoboken, NJ, 1 edition, 3 1990.  ISBN 978-0-471-87876-6.  doi:10.1002/9780470316801.
  * Patella and Ciaccia [2009] Marco Patella and Paolo Ciaccia.  Approximate similarity search: A multi-faceted problem.  _Journal of Discrete Algorithms_ , 7(1):36–48, 2009.  doi:10.1016/j.jda.2008.09.014.
  * Li et al. [2020] Pengfei Li, Hua Lu, Qian Zheng, Long Yang, and Gang Pan.  LISA: A Learned Index Structure for Spatial Data.  In David Maier, Rachel Pottinger, AnHai Doan, Wang-Chiew Tan, Abdussalam Alawini, and Hung Q. Ngo, editors, _Proceedings of the 2020 ACM International Conference on Management of Data, SIGMOD Conference 2020, online conference, Portland, OR, USA, June 14-19, 2020_ , SIGMOD ’20, page 2119–2133, New York, NY, USA, 2020. Association for Computing Machinery.  ISBN 9781450367356.  doi:10.1145/3318464.3389703.
  * Faloutsos and Lin [1995] Christos Faloutsos and King-Ip Lin.  FastMap: A Fast Algorithm for Indexing, Data-Mining and Visualization of Traditional and Multimedia Datasets.  In Michael J. Carey and Donovan A. Schneider, editors, _Proceedings of the 1995 ACM SIGMOD International Conference on Management of Data, San Jose, California, USA, May 22-25, 1995_ , pages 163–174. ACM Press, 1995.  doi:10.1145/223784.223812.
  * Duda et al. [2001] Richard O. Duda, Peter E. Hart, and David G. Stork.  _Pattern Classification_.  John Wiley & Sons, New York, NY, 2 edition, 2001.
  * Indyk and Motwani [1998] Piotr Indyk and Rajeev Motwani.  Approximate nearest neighbors: towards removing the curse of dimensionality.  In _Proceedings of the thirtieth annual ACM symposium on Theory of computing - STOC ’98_ , pages 604–613, Dallas, Texas, United States, 1998. ACM Press.  ISBN 978-0-89791-962-3.  doi:10.1145/276698.276876.
  * Beyer et al. [1999] Kevin Beyer, Jonathan Goldstein, Raghu Ramakrishnan, and Uri Shaft.  When Is “Nearest Neighbor” Meaningful?  In Catriel Beeri and Peter Buneman, editors, _Database Theory — ICDT’99_ , pages 217–235, Berlin, Heidelberg, 1999. Springer Berlin Heidelberg.  ISBN 978-3-540-49257-3.
  * Celebi et al. [2013] M. Emre Celebi, Hassan A. Kingravi, and Patricio A. Vela.  A comparative study of efficient initialization methods for the k-means clustering algorithm.  _Expert Systems with Applications_ , 40(1):200–210, January 2013.  ISSN 09574174.  doi:10.1016/j.eswa.2012.07.021.  URL http://arxiv.org/abs/1209.1960.
  * Aggarwal et al. [2001] Charu C. Aggarwal, Alexander Hinneburg, and Daniel A. Keim.  On the surprising behavior of distance metrics in high dimensional space.  In Gerhard Goos, Juris Hartmanis, Jan van Leeuwen, Jan Van den Bussche, and Victor Vianu, editors, _Database Theory — ICDT 2001_ , volume 1973, pages 420–434. Springer Berlin Heidelberg, Berlin, Heidelberg, 2001.  ISBN 978-3-540-41456-8.  doi:10.1007/3-540-44503-X_27.  Series Title: Lecture Notes in Computer Science.
  * Shvydun [2023] Sergey Shvydun.  Models of similarity in complex networks.  _PeerJ Computer Science_ , 9:e1371, May 2023.  ISSN 2376-5992.  doi:10.7717/peerj-cs.1371.  URL https://peerj.com/articles/cs-1371.
  * Ehsani and Drabløs [2020] Rezvan Ehsani and Finn Drabløs.  Robust Distance Measures for k NN Classification of Cancer Data.  _Cancer Informatics_ , 19:117693512096554, January 2020.  ISSN 1176-9351, 1176-9351.  doi:10.1177/1176935120965542.  URL http://journals.sagepub.com/doi/10.1177/1176935120965542.
  * Aumüller et al. [2018] Martin Aumüller, Erik Bernhardsson, and Alexander John Faithfull.  Ann-benchmarks: A benchmarking tool for approximate nearest neighbor algorithms.  _CoRR_ , abs/1807.05614, 2018.  URL http://arxiv.org/abs/1807.05614.
  * Duin and Pekalska [2012] Robert P.W. Duin and Elzbieta Pekalska.  The dissimilarity space: Bridging structural and statistical pattern recognition.  _Pattern Recognition Letters_ , 33(7):826–832, 2012.  ISSN 0167-8655.  doi:https://doi.org/10.1016/j.patrec.2011.04.019.  Special Issue on Awards from ICPR 2010.
  * Arthur and Vassilvitskii [2007] David Arthur and Sergei Vassilvitskii.  k-means++: the advantages of careful seeding.  In Nikhil Bansal, Kirk Pruhs, and Clifford Stein, editors, _Proceedings of the Eighteenth Annual ACM-SIAM Symposium on Discrete Algorithms, SODA 2007, New Orleans, Louisiana, USA, January 7-9, 2007_ , pages 1027–1035. SIAM, 2007.  URL https://dl.acm.org/doi/10.5555/1283383.1283494.
  * Tan et al. [2019] Pang-Ning Tan, Michael Steinbach, Anuj Karpatne, and Vipin Kumar.  _Introduction to data mining_.  Pearson, New York, NY, 2 edition, 2019.  ISBN 978-0-13-312890-1.
  * Bock [2008] Hans-Hermann Bock.  Origins and extensions of the k-means algorithm in cluster analysis.  _Journal Electronique d’Histoire des Probabilités et de la Statistique_ , 4(2), December 2008.  URL https://www.jehps.net/Decembre2008/Bock.pdf.
  * Maranzana [1964] F. E. Maranzana.  On the Location of Supply Points to Minimize Transport Costs.  _OR_ , 15(3):261, September 1964.  ISSN 14732858.  doi:10.2307/3007214.  URL https://www.jstor.org/stable/3007214?origin=crossref.
  * Diday [1971] E. Diday.  Une nouvelle méthode en classification automatique et reconnaissance des formes la méthode des nuées dynamiques.  _Revue de Statistique Appliquée_ , 19(2):19–33, 1971.  URL http://www.numdam.org/item/RSA_1971__19_2_19_0/.
  * Diday [1973] Edwin Diday.  The dynamic clusters method in nonhierarchical clustering.  _International Journal of Computer & Information Sciences_, 2:61–88, 03 1973.  doi:10.1007/BF00987153.
  * Diday et al. [1974a] Edwin Diday, A. Schroeder, and Y. Ok.  The dynamic clusters method in pattern recognition.  pages 691–697, 01 1974a.
  * Bock [1970] Hans-Hermann Bock.  _Automatische Klassifikation_ , pages 36–80.  01 1970.  ISBN 978-3-540-04962-3.  doi:10.1007/978-3-642-88253-1_10.
  * Diday and Govaert [1977] Edwin Diday and Gerard Govaert.  Automatic classification with adaptive intervals.  _RAIRO Inf Comput Sci_ , 11:329–349, 01 1977.
  * Kaufmann and Rousseeuw [1987] Leonard Kaufmann and Peter Rousseeuw.  Clustering by means of medoids.  _Data Analysis based on the L1-Norm and Related Methods_ , pages 405–416, 01 1987.
  * Schubert and Rousseeuw [2019] Erich Schubert and Peter J. Rousseeuw.  Faster k-Medoids Clustering: Improving the PAM, CLARA, and CLARANS Algorithms.  In _Similarity Search and Applications_ , pages 171–187, Switzerland, 2019. Springer Cham.  ISBN 978-3-030-32047-8.  doi:10.1007/978-3-030-32047-8_16.
  * Schubert and Zimek [2019] Erich Schubert and Arthur Zimek.  Elki: A large open-source library for data analysis - elki release 0.7.5 "heidelberg", 2019.
  * Schubert and Lenssen [2022] Erich Schubert and Lars Lenssen.  Fast k-medoids Clustering in Rust and Python.  _Journal of Open Source Software_ , 7(75):4183, July 2022.  ISSN 2475-9066.  doi:10.21105/joss.04183.  URL https://joss.theoj.org/papers/10.21105/joss.04183.
  * Bentley [1975] Jon Louis Bentley.  Multidimensional binary search trees used for associative searching.  _Commun. ACM_ , 18(9):509–517, September 1975.  ISSN 0001-0782.  doi:10.1145/361002.361007.  URL https://doi.org/10.1145/361002.361007.
  * Dong et al. [2011] Wei Dong, Charikar Moses, and Kai Li.  Efficient k-nearest neighbor graph construction for generic similarity measures.  In _Proceedings of the 20th International Conference on World Wide Web_ , WWW ’11, page 577–586, New York, NY, USA, 2011. Association for Computing Machinery.  ISBN 9781450306324.  doi:10.1145/1963405.1963487.
  * Andoni and Indyk [2006] Alexandr Andoni and Piotr Indyk.  Locality-sensitive hashing scheme based on p-stable distributions.  In _Proceedings of the eighteenth annual ACM-SIAM symposium on Discrete algorithms_ , pages 857–865. SIAM, 2006.  doi:10.1145/997817.997857.
  * Jégou et al. [2011] Hervé Jégou, Matthijs Douze, Cordelia Schmid, and Patrick Pérez.  Product quantization for nearest neighbor search.  _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , 33(1):117–128, 2011.  doi:10.1109/TPAMI.2010.57.
  * Aumüller et al. [2020] Martin Aumüller, Erik Bernhardsson, and Alexander Faithfull.  ANN-Benchmarks: A benchmarking tool for approximate nearest neighbor algorithms.  _Information Systems_ , 87:101374, 2020.  ISSN 0306-4379.  doi:https://doi.org/10.1016/j.is.2019.02.006.
  * Johnson et al. [1986] William Johnson, Joram Lindenstrauss, and Gideon Schechtman.  Extensions of lipschitz maps into banach spaces.  _Israel Journal of Mathematics_ , 54:129–138, 06 1986.  doi:10.1007/BF02764938.