1. Introduction

2. Key Euclidean Distance Properties for Computational Optimization

3. Key Approaches for Optimizing Euclidean Distance Computation

3.1. Algorithm-Specific Approaches

Algorithm-specific approaches for optimizing Euclidean distance computation focus on modifying or designing algorithms to reduce the number of required distance calculations by leveraging the mathematical properties of Euclidean distance, thereby minimizing computational overhead while preserving accuracy in distance measurements [33].

Squared Euclidean Distance

One of the primary methods for reducing the computational load of Euclidean distance calculations involves algorithmic optimizations that minimize the number of required operations. A straightforward optimization is to eliminate the computationally expensive square root operation when only relative distances between points are needed, as in K-nearest neighbors (KNN) or K-means clustering. Since the square root is a monotonically increasing function, the ordering of distances remains unchanged, allowing the use of squared Euclidean distance instead [1]. This simplification is particularly beneficial for large datasets, as square root operations are more computationally intensive than basic arithmetic. Additionally, using squared Euclidean distance not only enhances computational efficiency in the K-means algorithm but also improves its mathematical tractability [31].

For instance, in clustering techniques, the Euclidean distance is integral to solving the Minimum Sum-of-Squares Clustering (MSSC) problem [22], where the objective is to partition a set of data points $X={\left\{{\mathbf{x}}_i\right\}}_{i=1}^m\subset {\mathbb{R}}^n$ into K clusters ${\left\{{\mathcal{C}}_k\right\}}_{k=1}^K$ such that the sum of the squared Euclidean distances between each point and the centroid of its assigned cluster is minimized, that is by minimizing the total within-cluster variance:

$$\underset{{\left\{{\mathcal{C}}_k\right\}}_{k=1}^K}{min}\sum _{k=1}^K\sum _{{\mathbf{x}}_i\in {\mathcal{C}}_k}d{({\mathbf{x}}_i,{\mathit{\mu }}_k)}^2,$$

(2)

where $d{({\mathbf{x}}_i,{\mathit{\mu }}_k)}^2$ denotes the squared Euclidean distance between the data point ${\mathbf{x}}_i$ and the centroid ${\mathit{\mu }}_k$ of cluster ${\mathcal{C}}_k$. From a computational perspective, the squared Euclidean distance $d{({\mathbf{x}}_i,{\mathit{\mu }}_k)}^2$ is computed by simply eliminating the computationally expensive square root operation from the standard Euclidean distance and can be considered as one of the possible optimization approaches:

$$d{({\mathbf{x}}_i,{\mathit{\mu }}_k)}^2={\parallel {\mathbf{x}}_i-{\mathit{\mu }}_k\parallel }^2=\sum _{j=1}^n{(x_{i,j}-{\mathit{\mu }}_{k,j})}^2.$$

(3)

Given that the square root function is monotonically increasing, minimizing the squared Euclidean distance $d{({\mathbf{x}}_i,{\mathit{\mu }}_k)}^2$ is mathematically equivalent to minimizing the Euclidean distance $d({\mathbf{x}}_i,{\mathit{\mu }}_k)$:

$$\underset{{\mathbf{x}}_i}{min}d{({\mathbf{x}}_i,{\mathit{\mu }}_k)}^2\iff \underset{{\mathbf{x}}_i}{min}d({\mathbf{x}}_i,{\mathit{\mu }}_k).$$

This equivalence ensures that the cluster assignments obtained by minimizing the squared distance are identical to those obtained by minimizing the unsquared distance. This elimination simplifies computations, especially for large datasets, as square root operations are more computationally intensive than basic arithmetic operations.

The squared Euclidean distance $d{({\mathbf{x}}_i,{\mathit{\mu }}_k)}^2$ is a quadratic function, which is convex with respect to each data point ${\mathbf{x}}_i$ individually when the centroid ${\mathit{\mu }}_k$ is fixed. However, the overall objective function (2) in the MSSC algorithms is not convex with respect to both the cluster assignments and centroids simultaneously. Therefore, the K-means algorithm can converge to a local minimum, which is not necessarily the global minimum [34]. To find solutions closer to the global minimum, more sophisticated MSSC algorithms are typically employed [3].

Despite the non-convexity of the problem, the differentiability of the squared Euclidean distance allows for the straightforward derivation of update rules within the K-means algorithm. Specifically, the centroid ${\mathit{\mu }}_k$ of a cluster ${\mathcal{C}}_k$ can be straightforward and efficiently computed as the mean of the data points within that cluster:

$${\mathit{\mu }}_k={\displaystyle \frac{1}{|{\mathcal{C}}_k|}}\sum _{{\mathbf{x}}_i\in {\mathcal{C}}_k}{\mathbf{x}}_i.$$

(4)

The K-means algorithm iteratively assigns each data point ${\mathbf{x}}_i$ to the cluster with the nearest centroid ${\mathit{\mu }}_k$, and then updates each centroid ${\mathit{\mu }}_k$ by calculating the mean of all data points assigned to that cluster, i.e., ${\mathit{\mu }}_k={\displaystyle \frac{1}{|{\mathcal{C}}_k|}}{\sum }_{{\mathbf{x}}_i\in {\mathcal{C}}_k}{\mathbf{x}}_i$, iteratively minimizing the within-cluster variance until convergence.

A remarkable feature of the K-means algorithm is its ability to naturally and iteratively minimize the sum of squared Euclidean distances, ${\sum }_{k=1}^K{\sum }_{{\mathbf{x}}_i\in {\mathcal{C}}_k}{\parallel {\mathbf{x}}_i-{\mathit{\mu }}_k\parallel }^2$, making it particularly well-suited for large-scale clustering tasks. This natural minimization occurs through the algorithm’s intrinsic process of reducing the sum of squared distances between data points and their respective cluster centroids by alternating between assignment and update steps, ultimately leading to convergence without the need for external adjustments [34].

By utilizing the squared Euclidean distance (3), the K-means algorithm achieves computational efficiency and mathematical tractability, even though it may converge to a local minimum [31]. This approach enables the algorithm to find meaningful cluster configurations in a computationally feasible manner, though it may not always guarantee the global optimum of the clustering objective. However, this disadvantage can be mitigated by using simple techniques like Multi-start K-means [34] clustering, which helps to obtain near-optimal solutions.

Lower Bound Techniques

The block vector approximation [2] allows the K-means algorithm to avoid unnecessary full distance calculations by providing a lower bound on the Euclidean distance between a point and a cluster center. This lower bound enables the algorithm to make early decisions about whether a cluster center can possibly be the closest center, thereby eliminating the need for a complete calculation of computationally expensive Euclidean distance in many cases.

Initially, to form the block vector representation, the given a high-dimensional data point $\mathbf{x}\in {\mathbb{R}}^d$, it is subdivided into smaller sub-vectors, or blocks, of fixed sizes along the dimensions:

$$\mathbf{x}=\left({\mathbf{x}}^{\left(1\right)},{\mathbf{x}}^{\left(2\right)},\cdots ,{\mathbf{x}}^{\left(b\right)}\right),$$

where each ${\mathbf{x}}^{\left(i\right)}$ represents a subset of consecutive components of $\mathbf{x}$. For each block, the p-norm of the sub-vector is computed, resulting in a block vector ${\mathbf{x}}^B$:

$${\mathbf{x}}^B=\left(\parallel {\mathbf{x}}^{\left(1\right)}{\parallel }_p,\parallel {\mathbf{x}}^{\left(2\right)}{\parallel }_p,\cdots ,{\parallel {\mathbf{x}}^{\left(b\right)}\parallel }_p\right).$$

This process reduces the dimensionality of the original vector while retaining essential information about its structure. During the K-means algorithm, each point $\mathbf{x}$ is assigned to the closest cluster center. For a point $\mathbf{x}$ with a known closest center ${\mathbf{c}}_{\mathrm{current}}$, the lower bound on the distance to any other center ${\mathbf{c}}^{\prime }$ is computed using the block vectors. By using Hölder’s inequality, the inner product $\langle \mathbf{x},{\mathbf{c}}^{\prime }\rangle$ is approximated by the inner product of the block vectors $\langle {\mathbf{x}}^B,{\mathbf{c}}^{\prime B}\rangle$. If the lower bound on the distance to ${\mathbf{c}}^{\prime }$ is greater than the actual distance to the current closest center ${\mathbf{c}}_{\mathrm{current}}$, i.e.,

$$d(\mathbf{x},{\mathbf{c}}_{\mathrm{current}})\le \sqrt{{\parallel \mathbf{x}\parallel }_2^2+{\parallel {\mathbf{c}}^{\prime }\parallel }_2^2-2\langle {\mathbf{x}}^B,{\mathbf{c}}^{\prime B}\rangle },$$

then ${\mathbf{c}}^{\prime }$ cannot be the closest center, and the algorithm skips the full Euclidean distance calculation for ${\mathbf{c}}^{\prime }$. Here the expression ${\parallel \mathbf{x}\parallel }_2^2$ refers to the squared Euclidean norm of the vector $\mathbf{x}$. This results in significant computational savings, especially when there are many centers.

Algorithm 1 illustrates the snippet of K-means pseudocode leveraging lower bounds from block vector approximations to accelerate clustering.

Algorithm 1: Snippet of K-means pseudocode with Block Vectors optimization

The core computational saving comes from the fact that block vectors are much shorter than the original vectors, and the dot product of block vectors is much cheaper than a full Euclidean distance computation. Since ${\parallel \mathbf{x}\parallel }_2^2$ can be precomputed for all x, they don’t need to be recomputed in every iteration. This mechanism works particularly well when there are a large number of clusters, or when the data is high-dimensional, as the savings in computation increase in such scenarios.

Triangle Inequality

The triangle inequality property [1], $d({\mathbf{x}}_i,{\mathbf{x}}_k)\le d({\mathbf{x}}_i,{\mathbf{x}}_j)+d({\mathbf{x}}_j,{\mathbf{x}}_k)$, can be used to prune unnecessary distance calculations. In certain cases, if a partial distance is already larger than a known threshold, further calculation can be skipped, as the result will not affect the outcome [27].

For example, Hamerly’s method [35] for accelerating the K-means algorithm leverages the triangle inequality to significantly reduce the number of distance calculations required during clustering. In the naive K-means algorithm, the distance $d{({\mathbf{x}}_i,{\mathit{\mu }}_j)}^2$ between each data point ${\mathbf{x}}_i$ and every cluster centroid ${\mathit{\mu }}_j$ is recalculated in each iteration, resulting in a computational complexity of $O(m\times k)$ per iteration, where m is the number of points and k is the number of clusters. Hamerly’s method optimizes this process by maintaining upper bounds $u_i$ and lower bounds $l_i$ on the distances between each point ${\mathbf{x}}_i$ and its closest and second closest centroids, respectively. The algorithm also calculates a threshold $s_j$ for each cluster, representing half the minimum distance between the centroid ${\mathit{\mu }}_j$ and any other centroid: $s_j={\displaystyle \frac{1}{2}}{min}_{l\ne j}d({\mathit{\mu }}_j,{\mathit{\mu }}_l)$. Using the Triangle Inequality, it can skip unnecessary distance calculations if the current bounds indicate that the point cannot be closer to another centroid than its current assignment. Mathematically, this is achieved by comparing the upper bound $u_i$ to the maximum of $l_i$ and $s_j$, and only recomputing distances when necessary. Algorithm 2 presents a pseudocode snippet that demonstrates how Hamerly’s method efficiently reduces distance computations in the K-means algorithm when recalculating distances between points and centroids. In the presented pseudocode, $c_i$ denotes the index of the centroid (or cluster) currently assigned to data point ${\mathbf{x}}_i$. This approach effectively reduces the number of Euclidean distance calculations, resulting in a more efficient algorithm that converges faster than the naive K-means while maintaining the same accuracy in clustering.

Algorithm 2: Snippet of K-means pseudocode with Triangle Inequality optimization

Figure 1 illustrates the positions of data points and centroids on a number line, highlighting the distances and thresholds used to apply the triangle inequality in skipping distance computations during K-means clustering. In this example, the data point ${\mathbf{x}}_1$ is closer to centroid ${\mathit{\mu }}_1$ than to centroid ${\mathit{\mu }}_2$, and the upper bound $u_1$ is less than the threshold $s_1$. This condition allows us to skip the computation of the distance $d({\mathbf{x}}_1,{\mathit{\mu }}_2)$.

Recursive Distance Updating

Another widely-used method in data analysis that heavily relies on Euclidean distance computations is agglomerative hierarchical clustering [28], which builds a hierarchy by repeatedly merging the closest clusters until all points form a single cluster. Starting with each point as its own cluster, the algorithm merges the pair with the smallest distance and recalculates the distances between the newly merged cluster and all remaining clusters. In its naive implementation, this requires recalculating distances from scratch after every merge, making it computationally expensive for large datasets, as the number of pairwise distance calculations increases quadratically.

Ward’s method, one of the most efficient techniques for agglomerative hierarchical clustering, is an example of how to minimize Euclidean distance recalculations when solving the Minimum Sum-of-Squares Clustering (MSSC) problem. Ward’s method efficiently minimizes within-cluster variance (2) after merging clusters. It benefits from a recursive update mechanism known as the Lance–Williams formula, which optimizes the process of updating distances between the newly merged cluster and the remaining clusters. Instead of recalculating distances from scratch, the Lance–Williams formula updates the distances based on previously known distances and the sizes of the clusters involved in the merge. The Lance–Williams formula is:

$$D(i,z)={\displaystyle \frac{(n_i+n_x)D(i,x)+(n_i+n_y)D(i,y)-n_{i}D(x,y)}{n_x+n_y+n_i}}$$

Here, $D(i,z)$ is the distance between cluster i and the newly formed cluster z (created by merging clusters x and y); $n_i$, $n_x$, and $n_y$ are the sizes of clusters i, x, and y, respectively; $D(i,x)$, $D(i,y)$, and $D(x,y)$ are the pairwise squared Euclidean distances between these clusters.

Algorithm 3 presents pseudocode demonstrating how the Lance–Williams formula is used to efficiently update distances during agglomerative hierarchical clustering.

Algorithm 3: Hierarchical clustering algorithm using Lance-Williams formula for recursive distance updating

Using the Lance–Williams formula and the nearest-neighbor chain algorithm [36], agglomerative hierarchical clustering achieves a time complexity of $O\left(n^2\right)$, a significant improvement over the naive $O\left(n^3\right)$ approach.

Advanced Initialization

One effective strategy to reduce the computational burden in distance-intensive algorithms is to employ better initialization techniques, particularly in algorithms that rely on iterative local or global search [34]. Many such algorithms, like K-means clustering, repeatedly compute distances as they iteratively refine their solutions. Poor initialization, such as random selection of initial centroids in K-means, can lead to significantly more iterations and, consequently, a higher number of distance calculations. By contrast, advanced initialization methods can result in fewer iterations and thus fewer distance computations.

In K-means, for example, the algorithm starts by randomly selecting K centroids and iteratively refines them by recalculating distances between data points and centroids in each iteration. With random initialization, the algorithm often converges slowly due to poor starting points, requiring many iterations to find optimal centroids. Each iteration involves $O\left(nK\right)$ distance computations, where n is the number of data points and K is the number of clusters. Therefore, reducing the number of iterations directly reduces the number of distance calculations.

Advanced initialization techniques, such as K-means++ [37], provide a more strategic selection of initial centroids, leading to faster convergence. K-means++ initializes the centroids by probabilistically selecting data points that are farther apart, which ensures a better spread of centroids across the dataset. This reduces the number of iterations required for convergence compared to random initialization, ultimately reducing the number of Euclidean distance calculations performed during the algorithm.

An even more sophisticated initialization technique involves using the output of Ward’s method as the initial centroids for K-means. Ward’s method is aligns well with the objective of K-means [36]. By first applying Ward’s method and then using its cluster centroids as the starting points for K-means, the algorithm often converges in far fewer iterations than with K-means++ or random initialization. This combination can significantly reduce the number of distance computations, as the initial centroids are already near optimal. Furthermore, alternating between different MSSC clustering algorithms, such as using Ward’s method followed by K-means (i.e., using Ward’s output as the input for K-means), can often improve the clustering outcome by further minimizing the objective function (2) compared to running each algorithm in isolation, as the strengths of each method complement one another in refining the cluster assignments [22].

Precomputing and Caching Distances

In certain scenarios, such as when working with static datasets, it is possible to precompute the pairwise Euclidean distances between all points and store them in a distance matrix, which is a key component in many machine learning and data analysis tasks [26]. This precomputation, while expensive in terms of memory, allows for instant retrieval of distances during the optimization process, thereby eliminating the need for repeated calculations. Caching strategies can also be employed when the dataset is too large to store all pairwise distances, ensuring that frequently accessed distances are quickly available.

For a set of m points $\mathbf{X}=\{{\mathbf{x}}_1,{\mathbf{x}}_2,\cdots ,{\mathbf{x}}_m\}$ in n-dimensional space, the precomputed Euclidean distance matrix $\mathbf{D}$ is an $m\times m$ symmetric matrix, leveraging the property $d({\mathbf{x}}_i,{\mathbf{x}}_j)=d({\mathbf{x}}_j,{\mathbf{x}}_i)$, where the entry $D_{ij}$ represents the Euclidean distance between points ${\mathbf{x}}_i$ and ${\mathbf{x}}_j$ [32]. Mathematically, $\mathbf{D}$ is defined as:

$$D_{ij}=\left\{\begin{array}{cc}d({\mathbf{x}}_i,{\mathbf{x}}_j),\hfill & \mathrm{if}i>j\hfill \\ D_{ji},\hfill & \mathrm{if}i<j\hfill \\ 0,\hfill & \mathrm{if}i=j\hfill \end{array}\right.$$

(5)

In this formulation, the symmetry of the matrix is explicitly used by calculating $D_{ij}$ only for $i>j$ and setting $D_{ji}=D_{ij}$ for $i<j$. This approach allows for the computation of only half of the distance matrix—the elements below the main diagonal (the lower triangle), significantly reducing the number of distance calculations from $O\left(m^2\right)$ to $O\left(m\right(m-1)/2)$ for a dataset with m points.

To save space by storing only the lower triangular part of a symmetric matrix $\mathbf{D}$, its condensed (or compact) representation $\mathbf{C}$ can be used. The condensed matrix $\mathbf{C}$ is a vector containing the ${\displaystyle \frac{m(m-1)}{2}}$ elements of the lower triangular part of $\mathbf{D}$, excluding the diagonal. The index k in the condensed matrix $\mathbf{C}$, corresponding to the element $D_{ij}$ in the full matrix $\mathbf{D}$, is given by:

$$k(i,j,m)=i\times m+j-{\displaystyle \frac{i(i+3)}{2}}-1$$

where $0\le j<i\le m-1$ (assuming 0-based indexing). The condensed matrix $\mathbf{C}$ thus contains a total of ${\displaystyle \frac{m(m-1)}{2}}$ elements. To reconstruct the indices $(i,j)$ in the full matrix $\mathbf{D}$ from the index k in the condensed matrix $\mathbf{C}$, the following formulas can be used:

$$i=⌈{\displaystyle \frac{1}{2}}\left(-\sqrt{-8k+4m^2-4m-7}+2m-1\right)-1⌉,$$

$$j=⌊k-im+{\displaystyle \frac{i^2+3i+2}{2}}⌋$$

where $\lceil ·\rceil$ denotes the ceiling function, and $⌊·⌋$ denotes the floor function, which truncates the decimal part to return an integer value.

The bottleneck of the K-means algorithm is the repeated distance recalculations in each iteration. This is necessary because centroids are updated as the mean values of the assigned points, which are not actual data points. The K-means algorithm can be accelerated by using a precomputed Euclidean distance matrix $\mathbf{D}$, where $D_{ij}=d({\mathbf{x}}_i,{\mathbf{x}}_j)$, to avoid redundant computations. An optimized alternative is the K-medoids algorithm, which improves efficiency by leveraging the precomputed distance matrix and reducing the impact of outliers by using actual data points (medoids) as cluster centers [38].

Unlike K-means, which computes distances to centroids, K-medoids assigns points to medoids based on precomputed distances from matrix $\mathbf{D}$. The objective of K-medoids is to minimize the total dissimilarity between points and their medoids, defined as:

$$\underset{{\left\{{\mathbf{m}}_k\right\}}_{k=1}^K}{min}\sum _{k=1}^K\sum _{{\mathbf{x}}_i\in {\mathcal{C}}_k}d({\mathbf{x}}_i,{\mathbf{m}}_k)$$

where ${\mathbf{m}}_k\in {\mathcal{C}}_k$ is the medoid, and $d({\mathbf{x}}_i,{\mathbf{m}}_k)$ is the precomputed Euclidean distance. The algorithm iteratively refines the medoids by swapping them with non-medoid points and recalculating the total cost until convergence, making it especially beneficial for large datasets where recalculating distances is expensive.

In large datasets, there is a trade-off between computation time and memory usage when calculating Euclidean distances. Precomputing all pairwise distances in a matrix $\mathbf{D}$ speeds up algorithms by avoiding redundant recalculations, but this approach demands substantial memory. Storing the full distance matrix for m points requires $O\left(m^2\right)$ memory, which becomes impractical for large datasets. For example, with $m=100,000$, it would need around 80 GB of memory. While manageable for small datasets, this approach is infeasible for larger ones. To reduce memory usage, a condensed matrix $\mathbf{C}$ can be employed, leveraging the symmetry $D_{ij}=D_{ji}$ and zero diagonal $D_{ii}=0$, storing only the necessary elements.

Spatial Data Structures and Approximate Methods

The Euclidean distance is central in the K-nearest neighbors (KNN) algorithm, where it is employed to quantify the proximity between data points [23]. Given a set of training data ${\left\{({\mathbf{x}}_i,y_i)\right\}}_{i=1}^m$ where ${\mathbf{x}}_i\in {\mathbb{R}}^n$ represents the feature vectors and $y_i\in \mathcal{C}$ denotes the corresponding class labels, the KNN algorithm classifies a new data point $\mathbf{x}$ by identifying the set ${\mathcal{N}}_K\left(\mathbf{x}\right)$ of the K nearest neighbors, defined as:

$${\mathcal{N}}_K\left(\mathbf{x}\right)=\left\{{\mathbf{x}}_j\mid {\mathbf{x}}_j\in arg\underset{{\mathbf{x}}_i\in {\mathbb{R}}^n}{min}d(\mathbf{x},{\mathbf{x}}_i),j=1,2,\cdots ,K\right\},$$

The algorithm then assigns the class label y to $\mathbf{x}$ based on the majority class among its nearest neighbors:

$$y=arg\underset{c\in \mathcal{C}}{max}\sum _{{\mathbf{x}}_j\in {\mathcal{N}}_K\left(\mathbf{x}\right)}\mathbb{I}(y_j=c),$$

where $\mathbb{I}(·)$ is the indicator function. Thus, the Euclidean distance serves as a crucial metric in determining the classification of data points based on the proximity of their neighbors within the feature space.

Nearest neighbor search is a common task in many machine learning algorithms [6]. In nearest neighbor search the goal is to identify the closest point in a dataset $X={\left\{{\mathbf{x}}_i\right\}}_{i=1}^m\subset {\mathbb{R}}^n$ to a query point $\mathbf{x}\in {\mathbb{R}}^n$. The naive approach involves computing the Euclidean distance from $\mathbf{x}$ to each point in the dataset, which has a time complexity of $O\left(mn\right)$. For large datasets and high-dimensional spaces, this approach can be prohibitively expensive.

Besides the use of precomputed pairwise distance matrices, one of the most effective strategies to accelerate nearest neighbor search is the use of spatial data structures such as KD-trees [29] and Ball Trees [30]. These structures enable a partitioning of the dataset into hierarchical regions, allowing the algorithm to discard large portions of the dataset that are unlikely to contain the nearest neighbor, thus reducing the number of distance computations.

Spatial data structures, such as KD-trees, exploit the Spatial Locality property of Euclidean distance by organizing data points in a way that nearby points are grouped together and distant points are separated. The Spatial Locality property implies that points closer in space have smaller Euclidean distances, while those farther apart have larger distances.

The KD-tree algorithm operates by recursively splitting the dataset into two halves along the median of one of the dimensions at each level of the tree [29]. Given a query point $\mathbf{x}$, the KD-tree can be traversed to rapidly exclude regions of the space that are farther from $\mathbf{x}$ than the currently identified nearest neighbor. While the KD-tree is efficient in low to moderate dimensions, its performance degrades in high-dimensional spaces due to the curse of dimensionality.

For example, when searching for the nearest neighbor to a query point $\mathbf{x}$, the KD-tree first identifies the region where $\mathbf{x}$ resides. It then compares points within that region before expanding the search to neighboring regions only if necessary. This process avoids the need to compute distances to points in far-away regions, as the Euclidean distance between $\mathbf{x}$ and points in those distant regions would exceed the distances already found in the current local region. This selective search process dramatically reduces the computational burden, especially in low- and moderate-dimensional spaces.

Figure 2 illustrates how a KD-tree partitions the space and reduces distance computations during nearest neighbor search. The query point $\mathbf{x}$ is located in a specific region, and the algorithm only needs to consider data points within that region and nearby regions that intersect the search radius. Pruned regions are shaded, showing how the KD-tree efficiently eliminates the need to compute distances to distant points. Building spatial data structures like KD-trees or Ball Trees can require significant memory, especially with large datasets or higher dimensions.

An alternative approach is the use of Approximate Nearest Neighbors (ANN) algorithms, which trade-off exactness for speed. Techniques such as Locality-Sensitive Hashing (LSH) [39] enable faster searches by hashing data points into buckets based on random projections, ensuring that nearby points in the original space are more likely to be hashed into the same bucket. Although ANN algorithms do not guarantee the identification of the exact nearest neighbor, they offer significant computational savings and often suffice in practical applications.

Early Exit Strategies

Further optimizations can be achieved by leveraging techniques like early exit strategies in algorithms that do not require exact distances for all points. For instance, the additivity property of Euclidean distance in subspaces allows it to be broken down into the sum of squared differences across individual dimensions. This allows for partial distance calculations to determine if a point can be excluded early [40]. If the partial sum of squared differences already exceeds a known threshold, further dimensions do not need to be considered. This is an effective optimization in high-dimensional spaces, where early termination of distance computation can save substantial computational effort.

Dimensionality Reduction Techniques

As datasets grow in size and dimensionality, traditional methods like Euclidean distance face increasing computational challenges. The curse of dimensionality describes how the volume of space increases exponentially with the number of dimensions, causing data points to become sparse and distances less meaningful [24]. This sparsity complicates tasks such as nearest neighbor searches, clustering, and classification, as computational costs escalate. In high-dimensional spaces, Euclidean distance often fails to effectively differentiate between nearby and distant points, prompting the need for alternative distance measures and mitigation techniques [41].

To overcome these challenges, dimensionality reduction techniques are among the best options. These techniques aim to find a mapping $f:{\mathbb{R}}^n\to {\mathbb{R}}^d$, where $d\ll n$, that preserves the essential geometric properties of the data while reducing the computational cost of distance calculations.

Principal Component Analysis (PCA) [42] is one of the most widely used linear dimensionality reduction techniques. PCA transforms the original data points into a set of orthogonal components, ranked by the amount of variance they capture from the data. By retaining only the top d components, where $d\ll n$, PCA effectively reduces the dimensionality of the data while preserving the most significant features, allowing for faster Euclidean distance computations [43].

Another popular method is the t-Distributed Stochastic Neighbor Embedding (t-SNE) [44], which is a non-linear technique designed to maintain the local structure of the data in a lower-dimensional space. While t-SNE is primarily used for visualization, it can also facilitate faster distance computations in the lower-dimensional embedding space, although with the trade-off of potentially altering global distances.

Random Projection [45], grounded in the Johnson-Lindenstrauss Lemma, offers a probabilistic approach to dimensionality reduction. It projects the original data $\mathbf{X}\subset {\mathbb{R}}^n$ onto a randomly chosen lower-dimensional subspace ${\mathbb{R}}^d$ using a random matrix $\mathbf{R}\in {\mathbb{R}}^{n\times d}$, where each entry $R_{ij}$ is independently drawn from a suitable distribution (e.g., Gaussian or sparse distribution):

$$\mathbf{Z}={\displaystyle \frac{1}{\sqrt{d}}}\mathbf{XR}.$$

This projection approximately preserves pairwise Euclidean distances with high probability:

$$(1-ϵ)\left|\right|{\mathbf{x}}_i-{\mathbf{x}}_j{\left|\right|}^2\le \left|\right|{\mathbf{z}}_i-{\mathbf{z}}_j{\left|\right|}^2\le (1+ϵ)\left|\right|{\mathbf{x}}_i-{\mathbf{x}}_j{\left|\right|}^2,$$

where $ϵ\in (0,1)$ is a small distortion parameter. Random Projection is particularly appealing due to its computational efficiency and simplicity, as it requires no prior knowledge of the data distribution and can be implemented with minimal overhead.

Dimensionality reduction techniques offer computational advantages but come with limitations. PCA, though effective for preserving variance, is a linear method and may miss complex, non-linear patterns. Its sensitivity to outliers can distort the data structure, and it may discard low-variance dimensions that are still important. t-SNE, useful for visualizing local structure, is computationally expensive and can distort global relationships, making it unsuitable for large datasets. Random Projection is efficient and easy to implement but introduces distortion in distances, potentially affecting accuracy, and its results can vary due to randomness. These techniques often involve trade-offs between efficiency, accuracy, and data-specific considerations.

Approximation via Clustering and Dimensional Culling

Instead of computing the exact distance matrix, one can use clustering techniques to approximate the distance matrix by grouping similar points together. For example, using cluster centroids as representatives for the points within each cluster, the distance matrix can be computed more efficiently by calculating distances between centroids rather than between individual points, which significantly reduces the number of required distance calculations [46].

In high-dimensional spaces, not all dimensions contribute equally to the distance between points. Identifying and ignoring dimensions with low variance or limited impact on distance can reduce the number of calculations [24].

3.2. Low-Level Optimizations

Low-level optimizations focus on improving the efficiency of Euclidean distance computations by optimizing the underlying code at a granular level. These optimizations are crucial when working with large-scale data or real-time systems where every computational cycle counts.

Approximation with Lower Precision

By applying scalar quantization to the individual dimensions of a vector, each component of the vector is approximated using lower precision. For scalar quantization, each component of the vector is mapped to a discrete set of values. This reduces the complexity of the distance calculation by performing operations on smaller, quantized values rather than the original high-precision numbers. For instance, instead of using floating-point arithmetic, the vector components are represented as integers, allowing faster calculations using integer arithmetic, which is typically more efficient in many hardware architectures [47].

Scalar quantization maps each component $x_i\in \mathbb{R}$ of an n-dimensional vector $\mathbf{x}=(x_1,x_2,\cdots ,x_n)$ to a quantized value $q\left(x_i\right)\in \mathcal{Q}$, where $\mathcal{Q}$ is a finite set of levels. The resulting quantized value $q\left(x_i\right)$ is defined as $q\left(x_i\right)=arg{min}_{q_j\in \mathcal{Q}}|x_i-q_j|$. In uniform quantization, the levels in $\mathcal{Q}$ are equally spaced. For L levels over the range $[a,b]$, the quantization levels are given by $q_j=a+j·\Delta$, where $\Delta ={\displaystyle \frac{b-a}{L}}$ and $j=0,1,\cdots ,L-1$. The quantized value for $x_i$ is $q\left(x_i\right)=a+⌊{\displaystyle \frac{x_i-a}{\Delta }}⌉·\Delta$. The quantization error is the difference between $x_i$ and its quantized value: $e\left(x_i\right)=x_i-q\left(x_i\right)$. Although this error accumulates in distance calculations, it allows for faster computations using lower precision.

Quantization lowers the precision of vector components, introducing some error in distance calculations. The coarser the quantization, the larger the potential error. This trade-off is acceptable in cases where exact precision is less important, such as approximate nearest neighbor searches or noisy data. It is especially beneficial in large-scale or real-time applications, where faster computation and reduced memory use outweigh the accuracy loss.

Loop Unrolling

Loop unrolling is another low-level optimization technique that can be applied to Euclidean distance calculations. By manually unrolling loops, the number of loop control instructions (such as incrementing indices and checking loop termination conditions) is reduced, allowing the CPU to perform more useful work per cycle [48].

For example, consider the following loop that computes the squared Euclidean distance:

This loop can be unrolled by manually expanding the loop body for a fixed number of iterations:

While this increases the code size, it reduces the overhead of loop control and increases instruction-level parallelism, leading to faster execution.

Machine Code Optimization

Another avenue for optimizing Euclidean distance computations involves low-level programming and the use of just-in-time (JIT) compilers. Implementing distance calculations in low-level languages like C or C++ allows developers to exploit the full potential of the underlying hardware, using techniques such as assembly-level optimizations, manual loop unrolling, function inlining, instruction reordering, and precise control over memory allocation and access patterns to generate more efficient code [49]. These optimizations can minimize cache misses and take advantage of specific CPU instructions that are optimized for floating-point arithmetic. By enabling these optimizations, further performance gains can be achieved without requiring manual intervention. In environments where floating-point computation is expensive, using fixed-point arithmetic can additionally reduce the computational load.

In addition to traditional low-level languages, Python developers can achieve similar performance gains by using Numba, a high-performance JIT compiler [50]. Numba automatically translates Python functions into optimized machine code at runtime using the LLVM compiler library. This approach allows Python programs to reach speeds comparable to those of C or FORTRAN without requiring developers to leave the Python ecosystem. By simply applying a Numba decorator to a Python function, the function is automatically compiled into highly optimized machine code, often resulting in substantial speedups for numerical algorithms. This enables the efficient execution of distance calculations while retaining the flexibility and ease of use that Python offers.

Vectorization

Modern computing architectures offer vectorized instructions that can significantly accelerate distance calculations. Vectorization enables performing the same operation on multiple data elements simultaneously, leveraging the computational power of modern CPUs. Vectorized implementations and the use of advanced linear algebra libraries take advantage of special floating-point hardware, such as SIMD (Single Instruction, Multiple Data) operations and vector registers. By vectorizing the Euclidean distance computation, it is possible to leverage SIMD instructions to compute multiple distances in parallel within a single CPU core, thus speeding up the process [51].

Libraries such as BLAS (Basic Linear Algebra Subprograms) and NumPy provide highly optimized routines for these vectorized operations, making them readily accessible for optimization. For instance, BLAS library implementations include a set of low-level routines for performing common linear algebra operations—like matrix multiplication, dot products, vector addition, scalar multiplication, norm computation, and sum reduction—which can be directly used to calculate Euclidean pairwise distance matrices.

NumPy is particularly advantageous for vectorizing Euclidean distance computations due to its ability to efficiently handle large arrays and matrices through optimized low-level implementations [52]. These vectorized operations are built on highly optimized C libraries, allowing NumPy to fully leverage modern CPU architectures. By utilizing vectorization, NumPy can perform element-wise arithmetic across entire arrays simultaneously, significantly reducing the time complexity compared to traditional for-loop-based approaches in Python. This efficiency is further enhanced by the SIMD capabilities of modern CPUs, which enable the same operation to be applied to multiple data points in parallel. Additionally, NumPy’s integration with highly optimized linear algebra libraries, such as BLAS and LAPACK, ensures that operations like matrix multiplication and dot products—crucial for distance calculations—are executed with minimal overhead.

For example, vectorization techniques can be applied to distance calculations by expressing them in terms of matrix operations. The matrix of squared Euclidean distances between two vectors, ${\mathbf{x}}_i$ and ${\mathbf{y}}_j$, is given by:

$$D_{ij}=d{({\mathbf{x}}_i,{\mathbf{y}}_j)}^2=\parallel {\mathbf{x}}_i-{\mathbf{y}}_j{\parallel }^2=\parallel {\mathbf{x}}_i{\parallel }^2+{\parallel {\mathbf{y}}_j\parallel }^2-2\langle {\mathbf{x}}_i,{\mathbf{y}}_j\rangle$$

Here, $\parallel {\mathbf{x}}_i{\parallel }^2$ represents the squared norm of ${\mathbf{x}}_i$, $\parallel {\mathbf{y}}_j{\parallel }^2$ is the squared norm of ${\mathbf{y}}_j$, and $\langle {\mathbf{x}}_i,{\mathbf{y}}_j\rangle$ denotes the dot product between ${\mathbf{x}}_i$ and ${\mathbf{y}}_j$. Instead of performing full matrix multiplication, the squared norms $\parallel {\mathbf{x}}_i{\parallel }^2$ and $\parallel {\mathbf{y}}_j{\parallel }^2$ are computed more efficiently using row-wise operations. This optimization avoids unnecessary calculations, as only the squared norms for each row are required, significantly reducing the computational complexity.

Listing 1 provides an example of Python function for computing pairwise squared Euclidean distances between two sets of points using vectorization capabilities of NumPy library.

Listing 1: Python code for computing pairwise squared Euclidean distances between two sets of points using vectorization capabilities of NumPy

The vectorization capabilities of NumPy, as demonstrated in the code at Listing 1 and at Figure 3, significantly accelerate the computation of pairwise squared Euclidean distances by avoiding explicit Python loops and leveraging highly optimized matrix operations.

In particular, the dot product between matrices $\mathbf{X}$ and $\mathbf{Y}$ is computed using the function np.dot(X, Y.T), which takes advantage of NumPy’s efficient linear algebra routines to perform the operation in a single step across all pairs of points. Similarly, the squared norms of the rows in $\mathbf{X}$ and $\mathbf{Y}$ are computed using SIMD element-wise multiplication operation ($X*X$) and sum reduction along the row dimension (axis=1), fully utilizing vectorized operations to handle entire arrays at once. The operation $\mathtt{XX}[:,\mathtt{np}.\mathtt{newaxis}]$ reshapes the 1D array $XX$ into a 2D column vector, allowing NumPy to automatically expand the dimensions of $XX$ and $YY$ through broadcasting. In NumPy broadcasting enabling efficient element-wise operations with the 2D dot product result without explicitly copying data.

This approach eliminates the need for slow, iterative computations typically found in for-loop-based implementations. Instead, all operations are executed in parallel for each row of the matrices using special vector registers within a single CPU core, significantly boosting performance through SIMD operations. The final distance matrix is assembled by combining squared norms and dot products in a single vectorized expression, further optimizing computational efficiency.

Parallelization

Parallelization involves distributing the workload across multiple cores or processors, enabling simultaneous execution of distance calculations [3]. This approach is particularly advantageous when dealing with large datasets or when Euclidean distance computation needs to be performed repeatedly, as in iterative algorithms like K-means clustering [22] or K-nearest neighbors (KNN) [53]. By breaking down the computation into smaller tasks that can be processed concurrently, parallelization reduces the overall computation time, making it feasible to handle more extensive datasets within a reasonable timeframe.

Parallelization is a powerful strategy for optimizing Euclidean distance calculations, particularly when dealing with large-scale data or computationally intensive algorithms. By effectively utilizing the available hardware resources—whether through multi-threading, multiprocessing, or distributed computing—parallelization can dramatically improve the performance and scalability of algorithms that rely on Euclidean distance measurements.

The main conceptual difference between CPU and GPU parallelism lies in their design: CPUs are optimized for sequential processing with a few powerful cores, focusing on complex tasks that require strong single-thread performance, whereas GPUs are designed for massive parallelism with thousands of smaller, simpler cores that excel at handling many simultaneous tasks, making them ideal for workloads like distance computations across large datasets [53].

Libraries and frameworks like OpenMP for CPUs, CUDA for GPUs, and Apache Spark for distributed computing facilitate parallel computation by providing interfaces to efficiently utilize multiple cores, threads, and clusters, enhancing the performance of Euclidean distance calculations on various hardware architectures [54].

When combined with vectorization and machine code optimization, parallelization offers even greater potential for accelerating computations by leveraging modern hardware capabilities such as SIMD instructions, multi-core processors, and GPUs. By applying vectorized operations within each parallel thread or process, computational efficiency can be further enhanced, leading to significant reductions in processing time. Listing 2 provides an example of Python code for computing pairwise squared Euclidean distances between two sets of points using Numba’s JIT compilation. This approach combines automatic machine code optimization, vectorization, and parallelism across multiple CPU cores to enhance performance.

Listing 2: Python code for computing pairwise squared Euclidean distances between two sets of points using Numba’s JIT compilation, which combines automatic machine code optimization, vectorization, and parallelism across multiple CPU cores for enhanced performance.

In Listing 2, Numba’s JIT (Just-In-Time) compilation is applied to further accelerate the squared Euclidean distance computation compared to Listing 1. The @njit(parallel=True) decorator compiles the Python code into optimized machine code at runtime, allowing for faster execution by utilizing the CPU’s hardware capabilities. Parallelism is achieved through the prange construct, which distributes the workload of iterating over the rows of matrices X and Y across multiple CPU cores. This ensures that each core computes a portion of the pairwise distances, significantly speeding up the overall computation compared to the single-core, vectorized approach in Listing 1.

This parallelism significantly accelerates computation for larger datasets, especially when multiple CPU cores are utilized, as the pairwise distance calculations can be processed concurrently. As a result, the parallelization approach is particularly advantageous for large-scale problems, offering greater performance gains than vectorization alone, particularly when dealing with high-dimensional or large datasets and using multiple CPU or GPU cores. As the problem scale and the number of cores increase, the speedup achieved by parallelization can far exceed that of non-parallel implementations. However, when applied to small datasets, parallelization can increase memory usage due to the overhead of managing multiple threads or processes, as well as potential data duplication. Furthermore, combining vectorization and parallelization can yield extra efficiency by optimizing both individual arithmetic operations and the distribution of workload across cores, resulting in even greater overall acceleration.

Hardware Acceleration

Leveraging specialized hardware, such as multiple core CPUs and Graphics Processing Units (GPUs), can further accelerate Euclidean distance calculations [53]. Modern GPUs consist of thousands of cores that can execute the same operation on multiple data points simultaneously, making them well-suited for tasks like Euclidean distance calculation where the same operation (distance computation) needs to be applied across many data points. GPUs, with their massively parallel architecture, are particularly beneficial for applications involving large datasets, where the parallel computing power of GPUs can dramatically reduce computation time.

Field-Programmable Gate Arrays (FPGAs) offer a different kind of hardware acceleration by allowing the user to create custom circuits that are optimized for specific tasks, such as Euclidean distance computation [55]. Unlike GPUs, which are fixed-function devices, FPGAs can be reprogrammed to implement custom logic that executes the distance computations in an optimized fashion. FPGAs are particularly advantageous in applications where power efficiency is critical, as they can provide significant speedup with lower power consumption compared to GPUs or CPUs.

Figure 4. Vectorized and parallelized computation of pairwise squared Euclidean distances using Numba’s JIT compilation and parallel processing.

Figure 4. Vectorized and parallelized computation of pairwise squared Euclidean distances using Numba’s JIT compilation and parallel processing.

4. Comparative Analysis of Optimization Approaches

5. Discussion

6. Conclusion

Abbreviations

References

1. Deza, M.M.; Deza, E. Encyclopedia of Distances, 4 ed.; Springer: Berlin, Heidelberg, 2016; p. 756. [Google Scholar] [CrossRef]
2. Bottesch, T.; Bühler, T.; Kächele, M. Speeding up k-means by approximating Euclidean distances via block vectors. International conference on machine learning; PMLR, 2016; pp. 2578–2586. [Google Scholar]
3. Mussabayev, R.; Mussabayev, R. Superior Parallel Big Data Clustering Through Competitive Stochastic Sample Size Optimization in Big-Means. Intelligent Information and Database Systems; Nguyen, N.T., Chbeir, R., Manolopoulos, Y., Fujita, H., Hong, T.P., Nguyen, L.M., Wojtkiewicz, K., Eds.; Springer Nature Singapore: Singapore, 2024; pp. 224–236. [Google Scholar] [CrossRef]
4. Liberti, L.; Lavor, C. Euclidean distance geometry; Springer, 2017; Volume 3. [Google Scholar]
5. Croom, F.H. Principles of topology; Courier Dover Publications, 2016. [Google Scholar]
6. Braga-Neto, U. Fundamentals of pattern recognition and machine learning; Springer, 2020. [Google Scholar]
7. Oyewole, G.J.; Thopil, G.A. Data clustering: application and trends. Artificial Intelligence Review 2023, 56, 6439–6475. [Google Scholar] [CrossRef] [PubMed]
8. Varoquaux, G.; Buitinck, L.; Louppe, G.; Grisel, O.; Pedregosa, F.; Mueller, A. Scikit-learn: Machine learning without learning the machinery. GetMobile: Mobile Computing and Communications 2015, 19, 29–33. [Google Scholar] [CrossRef]
9. Burger, W.; Burge, M.J. Digital image processing: An algorithmic introduction; Springer Nature, 2022. [Google Scholar]
10. Tolebi, G.; Dairbekov, N.S.; Kurmankhojayev, D.; Mussabayev, R. Reinforcement learning intersection controller. 2018 14th International Conference on Electronics Computer and Computation (ICECCO); IEEE, 2018; pp. 206–212. [Google Scholar]
11. Fischer, M.M.; Scholten, H.J.; Unwin, D. Geographic information systems, spatial data analysis and spatial modelling: an introduction. In Spatial analytical perspectives on GIS; Routledge, 2019; pp. 3–20. [Google Scholar]
12. Tang, Y.; Zhao, L.; Zhang, S.; Gong, C.; Li, G.; Yang, J. Integrating prediction and reconstruction for anomaly detection. Pattern Recognition Letters 2020, 129, 123–130. [Google Scholar] [CrossRef]
13. Eiselt, H.A.; Sandblom, C.L. Decision analysis, location models, and scheduling problems; Springer Science & Business Media, 2013. [Google Scholar]
14. Carter, C.R.; Rogers, D.S.; Choi, T.Y. Toward the theory of the supply chain. Journal of supply chain management 2015, 51, 89–97. [Google Scholar] [CrossRef]
15. Perfilyeva, A.; Bespalova, K.; Kuzovleva, Y.; Mussabayev, R.; et al. . Genetic diversity and origin of Kazakh Tobet Dogs. Scientific Reports 2024, 14, 23137. [Google Scholar] [CrossRef]
16. van den Belt, M.; Gilchrist, C.; Booth, T.J.; Chooi, Y.H.; Medema, M.H.; Alanjary, M. CAGECAT: The CompArative GEne Cluster Analysis Toolbox for rapid search and visualisation of homologous gene clusters. BMC bioinformatics 2023, 24, 181. [Google Scholar] [CrossRef]
17. Mussabayev, R. Colour-based object detection, inverse kinematics algorithms and pinhole camera model for controlling robotic arm movement system. 2015 Twelve International Conference on Electronics Computer and Computation (ICECCO); IEEE, 2015; pp. 1–9. [Google Scholar]
18. Mukhamediev, R.I.; Yakunin, K.; Aubakirov, M.; Assanov, I.; Kuchin, Y.; Symagulov, A.; Levashenko, V.; Zaitseva, E.; Sokolov, D.; Amirgaliyev, Y. Coverage path planning optimization of heterogeneous UAVs group for precision agriculture. IEEE Access 2023, 11, 5789–5803. [Google Scholar] [CrossRef]
19. Bojanowski, P.; Grave, E.; Joulin, A.; Mikolov, T. Enriching word vectors with subword information. Transactions of the association for computational linguistics 2017, 5, 135–146. [Google Scholar] [CrossRef]
20. Li, J.; Wu, L.; Hong, R.; Hou, J. Random walk based distributed representation learning and prediction on social networking services. Information Sciences 2021, 549, 328–346. [Google Scholar] [CrossRef]
21. Altman, N.; Krzywinski, M. The curse (s) of dimensionality. Nat Methods 2018, 15, 399–400. [Google Scholar] [CrossRef]
22. Mussabayev, R.; Mladenovic, N.; Jarboui, B.; Mussabayev, R. How to use K-means for big data clustering? Pattern Recognition 2023, 137, 109269. [Google Scholar] [CrossRef]
23. Uddin, S.; Haque, I.; Lu, H.; Moni, M.A.; Gide, E. Comparative performance analysis of K-nearest neighbour (KNN) algorithm and its different variants for disease prediction. Scientific Reports 2022, 12, 6256. [Google Scholar] [CrossRef]
24. Aggarwal, C.C.; Hinneburg, A.; Keim, D.A. On the surprising behavior of distance metrics in high dimensional space. Database theory—ICDT 2001: 8th international conference, 2001 proceedings 8. London, UK, 4-6 January 2001; Springer, 2001; pp. 420–434. [Google Scholar]
25. Maitrey, S.; Jha, C. MapReduce: simplified data analysis of big data. Procedia Computer Science 2015, 57, 563–571. [Google Scholar] [CrossRef]
26. Qi, Z.; Xiao, Y.; Shao, B.; Wang, H. Toward a distance oracle for billion-node graphs. Proceedings of the VLDB Endowment 2013, 7, 61–72. [Google Scholar] [CrossRef]
27. Elkan, C. Using the triangle inequality to accelerate k-means. In Proceedings of the 20th international conference on Machine Learning (ICML-03); 2003; pp. 147–153. [Google Scholar]
28. Contreras, P.; Murtagh, F. Hierarchical clustering; 2015; pp. 103–124. [Google Scholar] [CrossRef]
29. Bentley, J.L. Multidimensional binary search trees used for associative searching. Commun. ACM 1975, 18, 509–517. [Google Scholar] [CrossRef]
30. Omohundro, S.M. Five Balltree Construction Algorithms. Technical Report TR-89-063; International Computer Science Institute, 1989. [Google Scholar]
31. Bock, H.H. Origins and extensions of the k-means algorithm in cluster analysis. Electronic journal for history of probability and statistics 2008, 4, 1–18. [Google Scholar]
32. Ying, Y.; Li, P. Distance metric learning with eigenvalue optimization. The Journal of Machine Learning Research 2012, 13, 1–26. [Google Scholar]
33. Cormen, T.H.; Leiserson, C.E.; Rivest, R.L.; Stein, C. Introduction to algorithms; MIT press, 2022. [Google Scholar]
34. Fränti, P.; Sieranoja, S. How much can k-means be improved by using better initialization and repeats? Pattern Recognition 2019, 93, 95–112. [Google Scholar] [CrossRef]
35. Hamerly, G. Making k-means even faster. In Proceedings of the 2010 SIAM International Conference on Data Mining (SDM); SIAM, 2010; pp. 130–140. [Google Scholar] [CrossRef]
36. Jeon, Y.; Yoon, S. Multi-Threaded Hierarchical Clustering by Parallel Nearest-Neighbor Chaining. IEEE Transactions on Parallel and Distributed Systems 2015, 26, 2534–2548. [Google Scholar] [CrossRef]
37. Arthur, D.; Vassilvitskii, S. k-means++: The advantages of careful seeding. In Proceedings of the eighteenth annual ACM-SIAM symposium on Discrete algorithms; SIAM, 2007; pp. 1027–1035. [Google Scholar]
38. Park, H.S.; Jun, C.H. A simple and fast algorithm for K-medoids clustering. Expert systems with applications 2009, 36, 3336–3341. [Google Scholar] [CrossRef]
39. Indyk, P.; Motwani, R. Approximate nearest neighbors: towards removing the curse of dimensionality. Proceedings of the Thirtieth Annual ACM Symposium on Theory of Computing; Association for Computing Machinery: New York, NY, USA, 1998; STOC’98; pp. 604–613. [Google Scholar] [CrossRef]
40. McNames, J. Rotated partial distance search for faster vector quantization encoding. IEEE Signal Processing Letters 2000, 7, 244–246. [Google Scholar] [CrossRef]
41. Beyer, K.; Goldstein, J.; Ramakrishnan, R.; Shaft, U. When is “nearest neighbor” meaningful? Database Theory—ICDT’99: 7th International Conference Jerusalem, Israel, 10-12 January 1999; Proceedings 7. Springer, 1999; pp. 217–235. [Google Scholar]
42. Jolliffe, I.T.; Cadima, J. Principal component analysis: a review and recent developments. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 2016, 374, 20150202. [Google Scholar] [CrossRef]
43. Gewers, F.L.; Ferreira, G.R.; Arruda, H.F.D.; Silva, F.N.; Comin, C.H.; Amancio, D.R.; Costa, L.d.F. Principal component analysis: A natural approach to data exploration. ACM Computing Surveys (CSUR) 2021, 54, 1–34. [Google Scholar] [CrossRef]
44. van der Maaten, L.; Hinton, G. Visualizing Data using t-SNE. Journal of Machine Learning Research 2008, 9, 2579–2605. [Google Scholar]
45. Bingham, E.; Mannila, H. Random projection in dimensionality reduction: applications to image and text data. Knowledge Discovery and Data Mining 2001. [Google Scholar] [CrossRef]
46. Borodin, A.; Ostrovsky, R.; Rabani, Y. Subquadratic approximation algorithms for clustering problems in high dimensional spaces. Machine Learning 2004, 56, 153–167. [Google Scholar] [CrossRef]
47. Rodriguez, A.; Segal, E.; Meiri, E.; Fomenko, E.; Kim, Y.J.; Shen, H.; Ziv, B. Lower numerical precision deep learning inference and training. Intel White Paper 2018, 3, 1–19. [Google Scholar]
48. Lee, S.; Gerstlauer, A. Data-dependent loop approximations for performance-quality driven high-level synthesis. IEEE Embedded Systems Letters 2017, 10, 18–21. [Google Scholar] [CrossRef]
49. Pikus, F.G. The Art of Writing Efficient Programs: An advanced programmer’s guide to efficient hardware utilization and compiler optimizations using C++ examples; Packt Publishing Ltd, 2021. [Google Scholar]
50. Lam, S.K.; Pitrou, A.; Seibert, S. Numba: A llvm-based python jit compiler. In Proceedings of the Second Workshop on the LLVM Compiler Infrastructure in HPC; 2015; pp. 1–6. [Google Scholar]
51. Patterson, D.A.; Hennessy, J.L. Computer Architecture: A Quantitative Approach, 5th ed.; The Morgan Kaufmann Series in Computer Architecture and Design, Elsevier Science & Technology; 2011. [Google Scholar]
52. Harris, C.R.; Millman, K.J.; Van Der Walt, S.J.; Gommers, R.; Virtanen, P.; Cournapeau, D.; Wieser, E.; Taylor, J.; Berg, S.; Smith, N.J.; others. Array programming with NumPy. Nature 2020, 585, 357–362. [Google Scholar] [CrossRef]
53. Masek, J.; Burget, R.; Karasek, J.; Uher, V.; Dutta, M.K. Multi-GPU implementation of k-nearest neighbor algorithm. 2015 38th International Conference on Telecommunications and Signal Processing (TSP); IEEE, 2015; pp. 764–767. [Google Scholar]
54. Dean, J.; Ghemawat, S. MapReduce: simplified data processing on large clusters. Commun. ACM 2008, 51, 107–113. [Google Scholar] [CrossRef]
55. Boikos, K.; Bouganis, C.S. A scalable fpga-based architecture for depth estimation in slam. International Symposium on Applied Reconfigurable Computing; Springer, 2019; pp. 181–196. [Google Scholar]
56. Gribel, D.; Vidal, T. HG-means: A scalable hybrid genetic algorithm for minimum sum-of-squares clustering. Pattern Recognition 2019, 88, 569–583. [Google Scholar] [CrossRef]