
%
\documentclass[11pt,reqno]{amsart}
\usepackage[cp1251]{inputenc}
\usepackage[english,russian]{babel}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{amsfonts}
\usepackage{graphicx}
\usepackage{xcolor}
\usepackage[colorlinks]{hyperref}
\usepackage{algorithm,algpseudocode}
\usepackage{tabularx}
\usepackage{array,multirow}
\newcommand{\argmin}{\mathop{\mathrm{argmin}}\limits}



\usepackage{titlesec}
\titlespacing{\section}{0pt}{2.5ex}{1.5ex}
\titlespacing{\subsection}{0pt}{1.5ex}{1ex}
\titlespacing{\subsubsection}{0pt}{1.5ex}{1ex}
\titleformat{\section}{\large\bfseries\centering}{\thesection}{1em}{}
\titleformat{\subsection}[runin]{\bfseries}{\thesubsection.}{0.5em}{}[.\mbox{\ }]
\titleformat{\subsubsection}[runin]{\bfseries}{\thesubsubsection.}{0.4em}{}[.\mbox{\ }]





%%%%%---------- ORCiD
\newcommand\orcidicon[1]{\href{https://orcid.org/#1}{\includegraphics[scale=0.02]{orcid.pdf}}}


\def\udcs{???.??} %Here the author places classificators of the paper according to Russian classification system
\def\mscs{??X??} %Here the author places  classificators according to the AMS classification list.
\setcounter{page}{144}

\newtheorem{lemma}{Lemma}
\newtheorem{theorem}{Theorem}
\newtheorem{definition}{Definition}
\newtheorem{corollary}{Corollary}
\newtheorem{proposition}{Proposition}



\hypersetup{linkcolor=blue,filecolor=black,urlcolor=blue, citecolor=blue}

\def\logo{\textcolor{blue}{{\bf\Huge S\raisebox{0.2ex}{\hspace{0.55ex}\raisebox{0.05ex}e\hspace{-1.65ex}$\bigcirc$}MR}}}

\def\semrtop
     {
  \vbox{
     \noindent\logo \hspace{-3mm}  \parbox{12cm}{     \begin{center}
     {\LARGE \textsf{ÑÈÁÈÐÑÊÈÅ ÝËÅÊÒÐÎÍÍÛÅ}} \\[2mm]
     {\LARGE \textsf{ÌÀÒÅÌÀÒÈ×ÅÑÊÈÅ ÈÇÂÅÑÒÈß}} \\[2mm]
     {\Large \textsf{Siberian Electronic Mathematical Reports}} \\[1mm]
     {\Large\tt{\textcolor{blue}{http://semr.math.nsc.ru}}}\\[1mm]
     \raisebox{1ex}{ISSN 1813-3304 }
     \vspace*{-3mm}
%     {\small Èíñòèòóò ìàòåìàòèêè èì. Ñ.Ë.Ñîáîëåâà ÑÎ ÐÀÍ}\\[-1mm]
%     {\small Sobolev Institute of Mathematics SB RAS}
     \end{center}
}
     \vspace{-1mm}
     \noindent
      \begin{tabular}{c}
     \hphantom{aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa!} \\
     \hline\hline
     \end{tabular}

     \vspace{3mm}
     {\flushleft\it Vol. 20, No. 2, pp. 144--145 (2023) \hspace{\fill}{\rm\small ÓÄÊ \udcs \hspace{-4mm}}}
     \newline  %Volume, pages in Russian, filled by Editorial board
     	{\rm\small \textcolor{blue}{https://doi.org/10.33048/semi.2023.20.???}}\hphantom{aaaaaaaaaaaaaaaaaaaaaaaa\;}{\rm\small MSC\ \ \mscs }
  }%\vbox
}




\begin{document}
\renewcommand{\refname}{References}
\renewcommand{\proofname}{Proof.}
\renewcommand{\figurename}{Fig.}




\thispagestyle{empty}

\title[NPU-accelerated clustering algorithms]{\large NPU-accelerated k-means and k-medoids clustering algorithms}


\author[I.L. VASILYEV]{{\bf I.L. Vasilyev}\protect\orcidicon{0000-0003-1862-4586}}
\author[T.V. GRUZDEVA]{{\bf T.V. Gruzdeva}\protect\orcidicon{0000-0002-4579-3927}}
\author[D.A. BOYARKIN]{{\bf D.A. Boyarkin}\protect\orcidicon{0000-0002-7048-2848}}
\author[A.V. USHAKOV]{{\bf A.V. Ushakov}\protect\orcidicon{0000-0001-8289-6266}} 



\address{Igor Leonidovich Vasilyev, 
\newline\hphantom{iii} Tatiana Vladimirovna Gruzdeva,
\newline\hphantom{iii} Anton Vladimirovich Ushakov,   % Êîíòàêòíûå äàííûå âñåõ àâòîðîâ è ìåñòî ðàáîòû óêàçûâàþòñÿ òîëüêî íà àíãëèéñêîì ÿçûêå
\newline\hphantom{iii} Matrosov Institute for System Dynamics and Control Theory of SB RAS, 
\newline\hphantom{iii} Lermontov str., 134,
\newline\hphantom{iii} 664033 Irkutsk, Russia}%
\email{\textcolor{blue}{\{vil,gruzdeva,aushakov\}@icc.ru}}%

\address{Denis Alexandrovich Boyarkin % Êîíòàêòíûå äàííûå âñåõ àâòîðîâ è ìåñòî ðàáîòû óêàçûâàþòñÿ òîëüêî íà àíãëèéñêîì ÿçûêå
\newline\hphantom{iii} Melentiev Energy Systems Research Institute of SB RAS, 
\newline\hphantom{iii} Lermontov str., 130,
\newline\hphantom{iii} 664033 Irkutsk, Russia}%
\email{\textcolor{blue}{denisboyarkin@isem.irk.ru}}%

\thanks{\sc Vasilyev,I.L., Gruzdeva, T.V.,   Boyarkin, D.A., Ushakov, A.V.,} 
\thanks{NPU-accelerated k-means and k-medoids clustering algorithms}
\thanks{\copyright \ 2025 Vasilyev I.L., Gruzdeva T.V.,   Boyarkin D.A., Ushakov A.V. }
\thanks{\rm The research is funded by the Ministry of Education and Science of the }
\thanks{\rm Russian
Federation (project 121041300065-9)}
\thanks{\it Received May, 31, 2025, Published December, 1, 2025}%


\semrtop \vspace{1cm}
\maketitle {\small

\vspace{-12pt}
\centerline{{\it Communicated by} {\sc }
%({\it filled by Editorial board})
}

\bigskip
\bigskip

\begin{quote}
\noindent{\bf Abstract:} Clustering is one of the basic tasks in unsupervised machine learning and data mining. The clustering process consists of partitioning data items from a dataset so that similar items are grouped into non-overlapping clusters. The so-called center-based clustering approach remains one of the most popular ones to formalize the clustering problem. In this case, partitioning data items into clusters performed by identifying the fixed number of so-called representatives (centers) of clusters. Then, clusters are formed by assigning each data item to the closest center. 

The $k$-means (minimum sum-of-squares) and $k$-medoids (mini\-mum sum-of-stars) problems are the most well-known center-based clustering problems and perhaps the most popular clustering mo\-dels at all. 
Despite their simplicity and popularity, $k$-medoids and $k$-means algorithms are both suitable for clustering relatively small datasets due to their high-computational complexity. The most time-consuming operation in these algorithms is the calculation of distances (measures of dissimilarity) between data items. 
To make the complex and consuming mathematical computations of such type, recently the Neural Processing Unit (NPU) was designed and optimized as a highly specialized processor to carry out AI challengers. 

We proposed and developed special variants of well-known $k$-medoids and $k$-means algorithms accelerated by NPU and parallel implementation.
The algorithms realized outperform not only the CPU but also the GPU in terms of performance, effectiveness and efficiency, which is proved by computational simulations.\medskip

\noindent{\bf Keywords:} clustering, neural processing unit, $k$-means, $k$-medoids, parallel computing, machine learning.
 \end{quote}
}

\bigskip

\section{Introduction}
Clustering is one of the basic tasks in unsupervised machine learning and data mining. The clustering process consists of partitioning data elements from a dataset so that similar elements are grouped into non-overlapping clusters. %Clustering has a wide range of applications, from marketing and customer segmentation \cite{TriBhaEsw2018} to image and speech recognition \cite{Japane2005,MaTuLuoW2021}, from bioinformatics \cite{AyyLin2015,FreyDueck2007} and computer vision \cite{GongPaw2015} to analysis of Ethereum smart contracts \cite{NorFiz2017}. 
The so-called center-based clustering approach remains one of the most popular ones to formalize the clustering problem. In this case, partitio\-ning data items into clusters is performed by identifying a fixed number of so-called cluster representatives or cluster centers. Then, clusters are obtained by assigning each data item to the closest (most similar) center. The number of clusters is equal to the number of centers, so it is supposed to be known in advance. Given a distance (or dissimilarity measure), any center-based clustering problem can be represented as a non-convex optimization problem, where the goal is to identify centers in such a way that the total sum of dissimilarities between data items and the closest centers is minimized. 

The $k$-means (minimum sum-of-squares) and $k$-medoids (minimum sum-of-stars) problems are the most well-known center-based clustering problems and perhaps the most popular clustering models at all. Note that they can be considered as facility location problems, e.g., $k$-medoids problem  is the well-known $p$-median facility location problem, and $k$-means model is a particular case of the multi-source Weber problem~\cite{VasUsh21,VasUshGru15}.

%The $k$-means and $k$-medoids clustering problems were first formulated as an optimization problem in~\cite{Vinod69} and still have been attracting much attention from different communities of researchers.  

Given a finite set of $m$ data items, represented as vectors $a^j\in I\!\! R^n$, $j\in J=\{1,\ldots,m\}$, the $k$-means problem is to find $k$ cluster centers (centroids) $x_c\in I\!\! R^n$, such that the total sum of squared Euclidean distances between data items and their closest centers is minimized:
\begin{equation}
\min_{C\subset I\!\! R^n}\Bigl\{\sum_{j=1}^{m}\min_{x_c\in C}\|a^j-x_c\|^2, \ |C|=k\Bigr\},\label{kmeans}
\end{equation} 
The problem is NP-hard even in the plane for arbitrary number of clusters $k$~\cite{MahNimVar12}. It is also NP-hard both in general dimension even for $k=2$~\cite{AloDesHanPop09} and when the dimension $n$ is a part of the input, whereas the number of clusters $k$ is not~\cite{DolKel11}.

%K-means and its modifications are the most popular algorithms for finding approximate solutions for the problem~\cite{kmeans}.

The $k$-medoids problem is very similar to the aforementioned $k$-means problem. It is to select $k$ medoids (or exemplars) $x_m^1,$ $x_m^2,$ $\ldots, x_m^k$ from within the finite set $\mathcal A$ so as  the overall sum of distances $d(\cdot,\cdot)$ (or other
measures of dissimilarity)  between data items and their closest medoids is minimized:
\begin{equation}
\min_{C\subset \mathcal A}\Bigl\{\sum_{j=1}^{m}\min_{x_m\in C}d(a^j,x_m), \ |C|=k\Bigr\}.\label{kmedoids}
\end{equation}
In contrast to the $k$-means problem, cluster representatives are not cluster means but so-called medoids selected among data items. The geometric (Euclidean and Manhattan) and graph versions of the $k$-medoids problem are known to be NP-hard \cite{KarivHakimi79,Megiddo79}.

Both clustering problems have widely been studied in the literature. One of the oldest and most popular clustering algorithms in practice, $k$-means (or Lloyd's algorithm~\cite{Lloyd82}), is a local search (alternative) heuristic to the minimum-sum-of-squares clustering problem. Besides $k$-means there are a number of other well-known algorithms, e.g. McQueen's algorithm and the Hartigan-Wong algorithm. In  spite of $k$-means, the latter is an exchange-type heuristic, i.e. it tries to improve the objective value by exchanging data items between clusters.  

The $k$-medoids problem is also known as the discrete $p$-median problem is one of the classical facility location problems. The most well-known clustering algorithm based on the $k$-medoids problem is PAM, which is a two stage heuristic including a greedy step to finds initial medoids followed by the best-improvement exchange heuristic over the SWAP neighborhood. The second stage is actually similar to the Hartigan-Wong algorithm. There is also an alternative (location-allocation) heuristic for the $k$-medoids problem, known as $k$-medoids or Cooper's algorithm, that follows the exactly same steps as $k$-means but finds medoids rather than means as cluster representative. Note that the main computational difference between k-means and $k$-medoids is the step of refining cluster centers. In case of $k$-means, for each cluster, it can be done in linear time resulting in the $\mathcal O(mnk)$ time complexity of a single iteration. For $k$-medoids, this step is not so simple and requires quadratic time in general case.      


%SWAP-based local search heuristics. K-medoids is a clustering algorithm, very similar to $k$-means, that is also known as Cooper's algorithm in case of the p-median problem.

Despite their simplicity and popularity, $k$-medoids and $k$-means are both suitable for clustering relatively small datasets due to their high-computatio\-nal complexity. For example, for MINIST8M test dataset, where $m = 8.1M$, $n = 784$, $k = 1000$, one iteration of $k$-means runs several hours~\cite{SheLiuTsa17}. However, real-worlds data may consists of much larger number of data items and possible clusters. Each data item may be defined by thousands of features, especially in applications involving clustering images or videos. Thus, it is extremely challenging to cluster such huge datasets with the conventional clustering algorithms.

There are several strategies of how make $k$-means and $k$-medoids applicable for large-scale data, e.g. random sampling (e.g. CLARA), binary hashing~\cite{SheLiuTsa17,GongPawYan15}, and fast procedures of distance calculations~\cite{DingZhaoShen15,YuChenChen20}. Note that as the most time consuming operations, finding the nearest cluster center and recompu\-ting cluster centers can be performed concurrently for each data item,% (cluster), 
there is a large strand of literature devoted the development of shared-memory and distributed parallel implementations of clustering algo\-rithms~\cite{UshVas21}. A huge number of computations per one iteration make GPUs especially suitable for speeding up the $k$-means and $k$-medoids algorithms making GPU-implemen\-tations the most effective and widespread. Effective GPU-accelerated imple\-mentations often requires the development of very efficient execution strategy. Thus, there are dozens of GPU-implementations for both clustering algo\-rithms that differ in what operations are performed on device or CPU, how to effectively implement data transfer, how to optimize basic algebraic operations, etc.~\cite{BayGhaHus18,BhiLeeMi15,CuoAngFarMar19,DaoZouMez19,GaoWuLiaMen24,KruKra20,LiFranPfa23,
LiZhaoChuLiu13,LutBreRab18,MaLim13,
MoWangXhang22,NelPal19,TayGow21}.%,WuDingZhaiLiu24}. 


The Neural Processing Unit (NPU) is highly specialized processor designed explicitly to carry out Artificial Intelligence (AI) challengers. It is optimized to make the complex and consuming mathematical computations
required in machine learning algorithms and artificial neural networks. Before the development of the NPU, all AI-related problems are solved by either the CPU or GPU (graphics processing unit) which are still effective and efficient, but not in AI area. NPUs excel in parallel processing capabilities, allowing to perform multiple operations simultaneously. This is crucial for solving problems that involve large datasets and complex computations. Second, focusing on energy-efficient designs, NPUs consume less power while delivering high performance, which make them ideally for edge devices. Furthermore, NPUs are tailored to accelerate specific AI tasks, e.g. deep learning, resulting in faster training and inference times.

In this paper we propose the first NPU-accelerated implementations of the well-known $k$-means and $k$-medoids algorithms. As was noted above the most time-consuming operation is the calculation of distances (measures of dissimilarity) between data items and cluster centers. Fortunately, these distances can be computed using special constructed matrices (tensors) and operations on them, which are performed very efficiently on the NPUs.  

% commentedto sent paper to zhangdong
%\cite{WuDingZhaiLiu24}
%\cite{MoWangXhang22}
%\cite{GaoWuLiaMen24}
%\cite{ZhangWangMo22}
%\cite{NelPal19}

%Yin Yam \cite{DingZhaoShen15} \cite{TayGow21}

%cuda k-means \cite{KruKra20}

%The rest of the paper is structured in the following manner. We give a brief description of the well-known $k$-means and $k$-medoids algorithms in Section 2. In Section 3, we present developed NPU-based implementations of these algorithms. In particular, we describe matrix representation and calculation on NPU the distance matrix, present the frameworks of the designed algorithms with Ascend NPU single operators' descriptions and features of Ascend memory management. The results of computational simulations are shown in Section 4.

%Besides being a basic data analysis tool, cluster analysis is a useful approach to image segmentation~\cite{ShaLiZhaYanSha19} or optimal location of service and distribution centers.





\section{Algorithms}
%The well-known $k$-means (Lloyd's) algorithm  remains one of the most popular clustering methods. It is an alternate heuristic similar to Cooper's algorithm~\cite{Coop64} proposed for the multi-source Weber problem. Its simplicity and relatively high efficiency for small size datasets made it popular in diverse applications. 
The $k$-means algorithm consists in two principal steps: assigning data items to the current closest (by the squared Euclidean distance) cluster centers followed by recompute the centers as means of data items assigned to the same clusters. The algorithm terminates when there are no changes in cluster assignments of a sufficiently large number of data items (see Alg.~\ref{alg:kmeans}). 
\begin{algorithm}[h]
\caption{$k$-means algorithm}\label{alg:kmeans}
\begin{algorithmic}[1]
\State Choose initial cluster centers $x_c^i\in\mathbb R^n$, $i\in I$ (often $x_c^i\in\{a^1\ldots,a^m\}$).
\While{clusters $C_i$ change}
\State Assign data items $a^j$, $j\in J$, to the closest cluster center $i^*$ by computing $i^* \gets \argmin_{i\in I}\sum\limits_{j=1}^m\|a^j-x_c^i\|^2$. Form clusters $C_i$ of data items assigned to the center $x_c^i$.
\State For each cluster $C_i$, $i\in I$, compute new $x_c^i$ as the center of gravity (mean) of all the data items from $C_i$: $x_c^i\gets \frac{1}{|C_i|}\sum\limits_{j:a^j\in C_i}a^j$.\label{alg:stepcenters}
\EndWhile
\end{algorithmic}
\end{algorithm}
%One can observe that an iteration of Alg.~\ref{alg:kmeans} requires $\mathcal O(kmn)$ time. The $k$-means algorithm is guaranteed to converge only to local optimal solutions. 
%Moreover, as was noted in~\cite{ArtVass07}, it may provide arbitrary bad partitions and has superpolynomial running time, e.g. for some datasets and specified initial cluster centers it requires $2^{\Omega(\sqrt{m})}$ iterations. Moreover, the running time remains superpolynomial with high probability even if the initial centroids are chosen uniformly at random~\cite{ArtVas06}. An upper bound on the running time for the case of the real line is $\mathcal O(m\Delta^2)$ where $\Delta$ is the ratio between the largest and the smallest distances between data items. 

The $k$-medoids algorithm \cite{ParkJun09} follows the same steps as $k$-means except for step~\ref{alg:stepcenters}. Instead of finding means, which can be done in linear time, $k$-meodis requires finding a medoid such that the distance between the data items in a cluster and their closet medoid is minimized. If on Step~\ref{alg:stepcenters} of Alg.~\ref{alg:kmeans}, for each cluster $C_i$, $i\in I$, we search for 
% compute new $y^i$ by search for 
$j_i \gets \argmin_{j:a^{j}\in C_i}\sum\limits_{\stackrel{j\prime:}{a^{j\prime}\in C_i}}\|a^j-a^{j\prime}\|^2$ and get $x_m^i\gets a^{j_i},$ then we obtain a $k$-medoids algorithm. Obviously, it can be done in quadratic time, which makes this step much more computationally expensive in comparison with $k$-means.

As $k$-medoids operates on the most centrally located objects in clusters, it is less sensitive to outliers in comparison with the $k$-means clustering.
%The algorithms differ mainly in Step 4 where new cluster's centers are calculated.
%If on Step 4 of Alg.~\ref{alg:kmeans}, for each cluster $C_i$, $i\in I$, 
%we search for  
% compute new $y^i$ by search for 
%$j_i \gets \argmin_{j:a^{j}\in C_i}\sum\limits_{\stackrel{j\prime:}{a^{j\prime}\in C_i}}\|a^j-a^{j\prime}\|^2$ and get $x_m^i\gets a^{j_i},$ then we obtain a $k$-medoids algorithm.
Moreover, $k$-medoids is flexible in the definition of dissimilarities between data items, e.g. dissimilarity measure may be not a proper distance metric. 

As $k$-means and $k$-medoids are local search heuristics, their performance are highly dependent on the strategy of selecting initial cluster centers.

\section{NPU-based implementation}
Here we present our implementations of the clustering algorithms. The main feature of NPUs that makes it highly effective in AI tasks is effective matrix computations. The idea of our NPU-based implementations of the clustering algorithms is to represent the location and allocation steps in the form of matrix-matrix computations. The latter allows us to achieve the full utilization of NPUs computing power.  %efficient work with matrix operations, them can be used to significantly accelerate clustering.
%In the algorithms for $k$-means and $k$-medoids clustering, the most time-consuming operation is calculation of distances (dissimilarities) between data items and centers of clusters that may be done by usage of matrix calculations.

\subsection{Matrix representation and calculations}
Let us first consider $k$-means. To form the matrix of distances between data items and cluster centers, $D\in I\!\! R^{m\times k}$, let us introduce the data matrix $A\in I\!\! R^{m\times n}$ where each row is the feature vector of a data item, and cluster center matrix $X\in I\!\! R^{k\times n}$ with rows corresponding to the coordinates of cluster centers:
$$
A=(a^1,a^2,\ldots,a^m)^\top,\;\;\;\; X=(x^1,x^2,\ldots,x^k)^\top.
$$
%	$$A=\left(
%	\begin{array}{cc}
%	a^1\\
%	a^2\\
%	\ldots\\
%	a^m\\
%	\end{array}
%	\right),\;\;\;\;
%	%$X\in I\!\!R^{k\times n}$ be a cluster centers' matrix, $
%	%\hspace{-0.8cm}
%	X=\left(
%	\begin{array}{cc}
%	x^1\\
%	x^2\\
%	\ldots\\
%	x^k\\
%	\end{array}
%	\right).
%	$$
Next, let us calculate squired Euclidean norm of each row and construct the following auxiliary matrices $A_{\|\cdot\|},X_{\|\cdot\|}\in I\!\! R^{m\times k} $:
%$A_d\in I\!\!R^{m\times k}$ be an item squared norms' matrix,
{%\small
\vspace{-5pt}
$$
%\hspace{-0.8cm}
A_{\|\cdot\|}=\left(
\underbrace{\begin{array}{ccc}
\parallel a^1\parallel^2&\ldots&\parallel a^1\parallel^2\\
\parallel a^2\parallel^2&\ldots&\parallel a^2\parallel^2\\
\ldots&\ldots&\ldots\\
\parallel a^m\parallel^2&\ldots&\parallel a^m\parallel^2\\
\end{array}}_{\Large k}
\right)\!\!,
%\hspace{-0.8cm}
X_{\|\cdot\|}=\left(\left.
\begin{array}{ccc}
\parallel x^1\parallel^2&\ldots&\parallel x^k\parallel^2\\
\parallel x^1\parallel^2&\ldots&\parallel x^k\parallel^2\\
\ldots&\ldots&\ldots\\
\parallel x^1\parallel^2&\ldots&\parallel x^k\parallel^2\\
\end{array}\right\} \!_m
\right)\!\!,
$$
}
Note that each row of the matrix $A_{\|\cdot\|}$ consists of $k$ copies of norms of the corresponding feature vector. Each column of the matrix $X_{\|\cdot\|}$ consists of $m$ copies of the corresponding cluster center. 

Thus, we can calculate the distance matrix $$
D=A_{\|\cdot\|}+X_{\|\cdot\|}-2AX^T,
$$
with elements $$
d_{ij}=\parallel a^j-x^i\parallel^2=\parallel a^j\parallel^2+\parallel x^i\parallel^2- 2\langle x^i,a^j\rangle,\;i=1,\ldots k,\;j=1,\ldots m.
$$

As for the $k$-medoids clustering, instead of one matrix, we compute the matrices $A_{i}\in I\!\! R^{m_i\times n} $ and $\widehat{A}_{\|\cdot\|}\in I\!\! R^{m_i\times m_i} $ for each cluster $i$ with $m_j$ data items:
\vspace{-5pt}$$A_i=\left(
\begin{array}{cc}
a^{i1}\\
a^{i2}\\
\ldots\\
a^{m_i}\\
\end{array}
\right),\;\;\;\;
%\hspace{-0.8cm}
\widehat{A}_{\|\cdot\|}=\left(
\underbrace{\begin{array}{ccc}
\parallel a^{i1}\parallel^2&\ldots&\parallel a^{i1}\parallel^2\\
\parallel a^{i2}\parallel^2&\ldots&\parallel a^{i2}\parallel^2\\
\ldots&\ldots&\ldots\\
\parallel a^{m_i}\parallel^2&\ldots&\parallel a^{m_i}\parallel^2\\
\end{array}}_{\Large m_i}
\right),
$$\\
Thus, for each cluster $i$, we can calculate 
the matrix of pairwise distances $D_i\in I\!\!R^{m_i\times m_i}:$ 
$$D_i=\widehat{A}_{\|\cdot\|}+\widehat{A}_{\|\cdot\|}^T-2{A}_i {A}_i^T$$ with elements 
$$
d_{ij}=\parallel a^i-a^j\parallel^2=\parallel a^i\parallel^2+\parallel a^j\parallel^2 - 2\langle a^i,a^j\rangle,\;\;\;i=1,\ldots m_i,\;j=1,\ldots m_i.
$$
Obviously, the pairwise distances between data items may be pre-computed and given to $k$-medoids as input. However, the large amount of memory required for storing these may be a serious bottleneck for the $k$-medoids algorithm, as their number is the square of the number of data items. %To avoid storing the whole distance matrix, one can apply a so-called $k$-nearest neighbors' strategy, so that for each data item $j$, one keeps in RAM only those $d_{ij}$, for which $i$ is among $k$-nearest neighbors of the item $j$. In this case, each column $j$ of $D$ is represented by the distances between data item $j$ and its $k$ nearest neighbors sorted in ascending order. 

\subsection{Huawei Ascend} 
The implementation of  $k$-means and $k$-medoids algorithms was carried out on the Ascend NPU devices using the Compute Architecture for Neural Networks (CANN). CANN is a device management framework that includes various APIs, particularly the AscendCL. AscendCL provides a collection of C language APIs to develop DNN (deep neural network) apps on the Ascend platform. AscendCL offers APIs for resource management, memory management, model loading and execution, operator loading and execution and other capabilities on the Ascend CANN platform with Ascend hardware. In a nutshell, AscendCL is a unified API framework used to invoke every resource. 

The core of Ascend's hardware computing capabilities is the so-called AI Processor, which efficiently performs matrix computations required in training and inference of neural networks, general computations, perform execution control, operations between scalars and vectors, etc. While these advantages are critical to performing DNN-related calculations, the strength of matrix calculations can be used in other algorithms that involve intensive matrix operations. In our implementation of the $k$-means and $k$-medoids algorithms, we focus on the attempt to exploit these advantages, developing matrix representations of location-allocation steps in the algorithms.

\subsection{Ascend memory management} 
The programs written using \linebreak AscendCL are executed on Devices. A Device refers to an Ascend AI Proces\-sor-powered hardware that connects to the Host over the PCIe interface, providing the Host with the NN computing capability. A Host is a computer on which one or more Devices are installed. The Device has its own installed RAM, which conceptually makes it similar to the GPU. In the case of using several Devices, their memory is independent of each other. Memory sharing between devices is not supported. The same can be said about the interaction between Device and Host.

Thus, performing any operations on Device means copying any data from the Host memory to the NPU Device for subsequent processing and copying the result back. This can be a costly procedure. On the other hand, limited installed memory on NPU Devices may require batch transfer of data for processing. These factors must be taken into account when developing algo\-rithms.

\subsection{NPU-accelerated $k$-means algorithm}\label{subs-means}
The implementation of $k$-means algorithm for NPU %using 
is shown in Figure \ref{fig1}, where the following notations are used:
\vspace{-6pt}
\begin{itemize}
   % \item $n$ -- dimension of feature vector;
    %\item $m$ -- number of data points;
    %\item $k$ -- number of clusters;
    \item $data$ -- ${m\times n}$ data points (host-side, C array);
    \item $max\_iter$ -- manually implemented iterations limit;
    \item $means$ -- ${k\times n}$ center matrix (host-side, C array);
    \item $A$ -- ${m\times n}$ data matrix (device-side, tensor);
    \item $X$ -- ${n\times k}$ center matrix (device-side, tensor);
    \item $P\!\_X$ -- ${n\times k}$ center matrix from the previous iteration (host-side, tensor);
    \item $PP\!\!\_X\!$ -- $\!{n\times k}$ center matrix from the iteration before previous (host-side, tensor);
    \item $A_{\|\cdot\|}$ -- ${m\times k}$ matrix with Euclidean norms of data items (host-side, tensor);
    \item $X_{\|\cdot\|}$ -- ${m\times k}$  matrix with Euclidean norms of cluster centers (host-side, tensor);
    \item $L$ -- ${m\times 1}$ point-to-cluster assignment labels (host-side, tensor).
\end{itemize}
\vspace{-6pt}
%The algorithm also uses other variables. Their description and dimension will be described in the algorithm itself. 
It is also worth mentioning that in terms of AscendCL, the concept of tensor is used. It is viewed as a multidimensional matrix. The $k$-means and $k$-medoids algorithms use ordinary two-dimensional matrices, therefore, %when describing the implementations of the algorithms below, 
a tensor should be understood here as a matrix.

    \begin{figure}[H]
    \centering
    \includegraphics[scale=0.7]{k-means alg.jpg}
    \caption{NPU-accelerated $k$-means algorithm scheme}
    \label{fig1}
\end{figure}

The algorithm uses AscendCL single operators and works as follows:
%\vspace{-6pt}
\begin{enumerate}
%\vspace{-3pt}
    \item Loading initial data matrix $data$ and initializing $means$ using the Lloyd algorithm;
    \item Transmitting the $data$ and $means$ to the memory of the Ascend NPU Device. The result is the tensors $A$ and $X$. To reduce the number of further transformations, $X$ is initialized in transposed form:\\
   \vspace{-5pt}  $$\begin{array}{l}
  data (C array) \rightarrow CreateAclTensor \rightarrow A (aclTensor) {m\times n};\\
means (C array) \rightarrow CreateAclTensor \rightarrow Y (aclTensor) {n\times k}.
  \end{array}$$
Memory is also reserved in the device for matrices: 
    $P\_X$, $PP\_X$ (centroid values at the previous iteration and the iteration before last), $X_{\|\cdot\|}$ (norms of the $X$ matrix), $L$ (pointer matrix). All further operations starting from this step are performed in the memory of the NPU Device and do not affect the computer's RAM, this saves resources on data transfer between Host and Device.
    \item The calculation of the tensor $A_{\|\cdot\|}$ is performed. % - the norm of matrix $A$. 
    Matrix $A_{\|\cdot\|}$ does not change during the calculation process, therefore it is calculated 
    %is performed 
    before starting  the iterative process only once. AscendCL provides a ready-made aclnnNorm operator to calculate the norm. However, for our problem, the result of its execution also needs to be squared, and the matrix dimension has been changed using the repetition operator $k$ times:
\vspace{-3pt}$$\begin{array}{l}
A \{m\times n\} \rightarrow aclnnNorm \rightarrow out\_temp \{M\times 1\} \rightarrow \\
    aclnnInplacePowTensorScalar(2) \rightarrow out\_temp \{M\times 1\} \rightarrow \\aclnnRepeatInterleaveIntWithDim(k) \rightarrow A_{\|\cdot\|} \{m\times k\}.\end{array}$$
    Changing the dimension by copying columns is necessary to perform further matrix operations.
 
    \item %Calculation of matrix $X_{\|\cdot\|}$ - norms of matrix $X$. 
    The calculation of the matrix $X_{\|\cdot\|}$ is identical to calculation of $A_{\|\cdot\|}$ (step~3), although it is performed for columns, not rows:
\vspace{-3pt}$$\begin{array}{l}
X \{n\times k\} \rightarrow  aclnnNorm \rightarrow  out\_temp \{1\times k\} \rightarrow  \\aclnnInplacePowTensorScalar(2) \rightarrow  out\_temp \{1\times k\} \rightarrow  \\aclnnRepeatInterleaveIntWithDim(m) \rightarrow  X_{\|\cdot\|} \{m\times k\}.\end{array}$$

    \item To calculating  the distance matrix $D$, AscendCL provides two single operators -- $aclnnAdd$ and $aclnnInplaceAddmm$. The first one is a regular summation, the second performs multiplication with a scalar and summation in one operation:
   \vspace{-3pt} 
   $$\begin{array}{l}
  A_{\|\cdot\|} \{m\times k\}, X_{\|\cdot\|} \{m\times k\} \rightarrow aclnnAdd \rightarrow  res\_temp \{m\times k\};\\
    res\_temp \{m\times k\}, A \{m\times n\}, X \{n\times k\} \rightarrow \\ aclnnInplaceAddmm(-2) \rightarrow  D \{m\times k\}.\end{array}$$
     
    \item Calculating the minimum distance from cluster's centers to each point and obtain a vector $L$ of pointers. To complete this step, one operator is enough:\\
  %\vspace{-3pt} 
   $D \{m\times k\} \rightarrow aclnnArgMinGetWorkspaceSize \rightarrow L \{M\times 1\}.$
 
    \item Saving the values of previous centroids:
  \vspace{-5pt} 
  $$\begin{array}{l} P\_X \{n\times k\} \rightarrow aclnnInplaceCopy \rightarrow PP\_X \{n\times k\};\\%[2pt]  
    X \{n\times k\} \rightarrow aclnnInplaceCopy \rightarrow P\_X \{n\times k\}.\end{array}$$
    
    \item Obtaining a binary matrix-mask reflecting the belonging of the points of matrix $A$ to the centroids $X$. To obtain it, the following operations are performed:
    \begin{enumerate}
        \item Transpose matrix $L$ and change its dimension using row repetiti\-on:
        \vspace{-5pt} 
  $$\begin{array}{l} L \{m\times 1\} \rightarrow aclnnPermute \rightarrow L^T \{1\times m\} \rightarrow \\aclnnRepeatInterleaveIntWithDim(k) \rightarrow L^T\!\_R \{k\times m\}.\end{array}$$
  \vspace{-5pt} 
        \item Creating a matrix of centroid numbers from array of cluster indexes:
\vspace{-5pt} 
$$\begin{array}{l}  Clusters(C array)\rightarrow CreateAclTensor \rightarrow \\clustersNum\_temp \{k\times 1\} \rightarrow \\
 aclnnRepeatInterleaveIntWithDim(m) \rightarrow \\clustersNumR\_temp \{k\times m\}.\end{array}$$
\vspace{-5pt} 
        \item Logical correlation of the contents of the pointer matrix and the matrix of centroid indexes allows us to obtain the required $mask$ matrix:
        \vspace{-3pt} 
$$\begin{array}{l}
L^T\!\_R \{k\times m\},\; clustersNumR\_temp \{k\times m\} \rightarrow \\aclnnEqTensor \rightarrow mask\_temp \{k\times m\}.\end{array}$$
        As a result we get a matrix of booleans. To use it in matrix operations below, conversion to float is needed:  
        \vspace{-5pt} 
$$\begin{array}{l}
         mask\_temp \{k\times m\} \rightarrow aclnnCast(float) \rightarrow \\mask\_temp \{k\times m\}.\end{array}$$
  \end{enumerate}    
  
  \item Calculating new centroids and updating the matrix $X$ need perfor\-ming the following steps:
  \begin{enumerate}
    
        \item Coordinate sums are calculated for all points included in the cluster:
\vspace{-5pt}      $$\begin{array}{l}  mask\_temp \{k\times m\}, \;A \{m\times n\} \rightarrow\\ aclnnMatmul \rightarrow  sum\_temp \{k\times n\}.\end{array}$$
        \item Calculation of the number of points included in each cluster:
\vspace{-5pt} $$\begin{array}{l}
mask\_temp \{k\times m\} \rightarrow aclnnReduceSum \rightarrow \\count\_temp \{k\times 1\} \rightarrow \\aclnnRepeatInterleaveIntWithDim(n) \rightarrow \\count\_R\_temp \{k\times n\}.\end{array}$$
        \item Calculation of new centroids:
\vspace{-5pt} $$\begin{array}{l}
sum\_temp \{k\times n\},\; count\_R\_temp \{k\times n\} \rightarrow\\ aclnnInplaceDiv \rightarrow X^T \{k\times n\} \rightarrow \\aclnnPermute \rightarrow  X \{n\times k\}.\end{array}$$

        The resulting matrix  $X$ may contain $NAN$ if one or more clusters are empty. For such clusters, the centroid is preserved from the previous iteration. In this case, search for $NAN$ is performed and  the coordinates of the corresponding centroids are replaced with the previous ones:
$$\begin{array}{l}
X \{n\times k\} \rightarrow aclnnIsFinite \rightarrow condition\_temp \{n\times k\};\\[4pt]
        condition\_temp \{n\times k\},\; X \{n\times k\},\;P\_X \{n\times k\} \rightarrow \\aclnnSWhere \rightarrow Y {N, K}.\end{array}$$
    \end{enumerate}
    \item The stopping criteria are checking. If the  matrix $X$ is equal to  $P\!\_X$ or $PP\!\_X$ matrices, then it means that the optimal centroids for the clusters have been found and the  matrix $X$ no longer changes, the cycle ends. Otherwise, repeat steps 4--9. Also, the terminate condition  is  to reach the maximum number of iterations.
    \item The matrix $X$ and labels $L$ are uploaded to the host and the result is output.
\end{enumerate}
\begin{figure}%[t!h]
    \centering
    \includegraphics[scale=0.65]{k-medoids alg.jpg} %[width=1\linewidth]
    \caption{$k$-medoids algorithm}
    \label{fig2}
\end{figure}
  \vspace{-6pt} 
  \subsection{$k$-medoids NPU algorithm}
The scheme  of the $k$-medoids algorithm implementation %in matrix representation for 
on NPU  is shown in Figure \ref{fig2}. The algorithm description uses the notations from   Subsection~\ref{subs-means} with the following additions and clarifications:
  \begin{itemize}
 %   \item $n$ - dimension of feature vector;
 %   \item $m$ - number of data points;
 %   \item $k$ - number of clusters;
 %   \item $data$ - ${m\times n}$ data points (host-side, C array);
 %   \item $max\_iter$ - manually implemented iterations limit;
    \item $medoids$ -- ${k\times 1}$ medoids' %matrix of center 
    indexes (host-side, C array);
   % \item $A$ - ${m\times n}$ data matrix (device-side, tensor);
    \item $Xmed$ -- ${k\times 1}$ center (medoids) indexes (device-side, tensor);
    \item $P\!\_Xmed$ -- ${k\times 1}$ medoids' %center matrix 
    indexes from the previous iteration (host-side, tensor);
    \item $PP\!\_Xmed$ -- ${k\times 1}$ medoids' %center matrix 
    indexes from the iteration before previo\-us (host-side, tensor);
%    \item $A_{\|\cdot\|}$ - ${m\times k}$ Euclidean norms of data matrix (host-side, tensor);
%    \item $X_{\|\cdot\|}$ - ${m\times k}$ Euclidean norms of center matrix (host-side, tensor);
 %   \item $L$ - ${m\times 1}$ point-to-cluster assignment labels (host-side, tensor).
\end{itemize}
%The algorithm also uses other variables. Their description and dimension will be described in the algorithm itself.



The algorithm uses AscendCL single operators and iterates as follows:
\begin{enumerate}
    \item Loading initial  $data$ matrix and initializing  $medoids$ % using the Lloyd algorithm 
    (pick $k$ random data items). The $medoids$ stores only indexes of vectors from $data$ accepted as medoids of clusters.
    \item Transmitting the $data$ and $medoids$ to the memory of the Ascend NPU Device. The result is the tensors $A$ and $Xmed$
      \vspace{-3pt} 
$$\begin{array}{l}
data (C array) \rightarrow CreateAclTensor \rightarrow A (aclTensor) \{m\times n\};\\%[4pt]
    medoids (C array) \rightarrow CreateAclTensor \rightarrow Xmed (aclTensor) \{k\times 1\}.
    \end{array}$$
    In the device, memory is also reserved  for matrices: 
    $P\!\_Xmed$,\linebreak $PP\!\_Xmed$ (medoids indexes at the previous iteration and the iterati\-on before last), 
    $X$(medoids), $X_{\|\cdot\|}$(norms of the medoids from $X$ matrix), $L$(pointer matrix).
    
    \item Calculating the tensor $A_{\|\cdot\|}$ (norms of the data items from matrix $A$). This step is completely identical to the step 3 of the $k$-means algorithm (see Subsection~\ref{subs-means}).%:

   % $A \{m\times n\} \rightarrow aclnnNorm \rightarrow out\_temp \{M\times 1\} \rightarrow 
  %  aclnnInplacePowTensorScalar(2) \rightarrow out\_temp \{M\times 1\} \rightarrow aclnnRepeatInterleaveIntWithDim(k) \rightarrow A_{\|\cdot\|} \{m\times k\}$;
    
   % Changing the dimension by copying columns is necessary to perform further matrix operations.
    \item Obtaining the matrix $X$ of the medoids:
\vspace{-6pt} $$\begin{array}{l}   
    A \{m\times n\}, \;Xmed \{k\times 1\} \rightarrow aclnnEmbedding \rightarrow \\X^T \{k\times n\} \rightarrow aclnnPermute \rightarrow X \{n\times k\}.
          \vspace{-4pt}\end{array}$$
    
    \item Calculating the matrix $X_{\|\cdot\|}$% - norms of matrix $X$. This step 
    is completely identical to the step 4 of the $k$-means algorithm (see Subsection~\ref{subs-means}).
    
  %  $X \{n\times k\} \rightarrow aclnnNorm \rightarrow out\_temp \{1\times k\} \rightarrow aclnnInplacePowTensorScalar(2) \rightarrow out\_temp \{1\times k\} \rightarrow aclnnRepeatInterleaveIntWithDim(m) \rightarrow X_{\|\cdot\|} \{m\times k\}$; 
   
    \item Calculating the distance matrix $D$ (see the AscendCL operators chain from the step 5 of the $k$-means algorithm).%:
    
   %$A_{\|\cdot\|} \{m\times k\}, X_{\|\cdot\|} \{m\times k\} \rightarrow aclnnAdd \rightarrow res\_temp \{m\times k\}$;
    
    %$res\_temp \{m\times k\}, A \{m\times n\}, X \{n\times k\} \rightarrow aclnnInplaceAddmm(-2) \rightarrow D \{m\times k\}$;
    
    \item Computing the minimum distance for each data item and obtaining a vector of pointers $L$. For further operating, $L$ is copied into  Host's RAM:
   \vspace{-4pt} 
   $$\begin{array}{l}  
 D \{m\times k\} \rightarrow aclnnArgMinGetWorkspaceSize \rightarrow \\L \{M\times 1\} \rightarrow  aclrtMemcpy \rightarrow clusters (C array).        \end{array}$$

    \item Saving the indexes of previous medoids:
 \vspace{-4pt}  $$\begin{array}{l} 
 P\!\_Xmed \{k\times 1\} \rightarrow aclnnInplaceCopy \rightarrow PP\!\_Xmed \{k\times 1\};\\%[4pt]
    Xmed \{k\times 1\} \rightarrow aclnnInplaceCopy \rightarrow P\!\_Xmed \{k\times 1\}.
    \end{array}$$ 
    \item The number of points $p$ in each cluster $cl$ is calculated, as well as the indices of these points in $data$ -- $P\!\_ind$. For this step, the vector $clusters$ is used. Operations are performed in Host's RAM. The resulting indices are copied into the memory of the Ascend NPU device:\\
$\;\;\;\;\;\;p\_ind (C array) \rightarrow CreateAclTensor \rightarrow P\!\_ind \{p\times 1\}.$
    \item A submatrix $A_P$ is formed by %, which includes 
    only the points under the indices from the set $P\!\_ind$:\\  
    $\;\;\;\;P\!\_ind \{p\times 1\},\; A \{m\times n\} \rightarrow aclnnEmbedding \rightarrow A_P \{p\times n\}.$\\[4pt]
    For the resulting matrix, the norms of data items are calculated:
 \vspace{-4pt} 
 $$\begin{array}{l} 
 A_P \{p\times n\} \rightarrow aclnnNorm \rightarrow out\_temp \{p\times 1\} \rightarrow \\aclnnInplacePowTensorScalar(2) \rightarrow out\_temp \{p\times 1\} \rightarrow \ \end{array}$$   
    For both matrices, a transpose operation is performed:
 \vspace{-4pt} $$\begin{array}{l} 
 A_P \{p\times n\} \rightarrow aclnnPermute \rightarrow A_P^T \{n\times p\};\\%[4pt]    
    A_{P\|\cdot\|} \{p\times p\} \rightarrow aclnnPermute \rightarrow A_{P\|\cdot\|}^T \{p\times p\}.
      \vspace{-6pt} \end{array}$$
    These matrices are needed to calculate the distance matrix in the next step.
    
   
    \item Calculating the distance submatrix $D_P$ by usage the matrices $A_{P\|\cdot\|} $ and $ A_{P\|\cdot\|}^T$. This step is identical to the step 5 of the $k$-means algorithm and the step 6 of this algorithm with precision to the notations.
    
   % $A_{P\|\cdot\|} \{p\times p\}, A_{P\|\cdot\|}^T \{p\times p\} \rightarrow aclnnAdd \rightarrow res\_temp \{p\times p\}$;
   %    $res\_temp \{p\times p\}, A_P \{p\times n\}, A_P^T \{n\times p\} \rightarrow aclnnInplaceAddmm(-2) \rightarrow D \{p\times p\}$;
   
    \item Finding the point with the smallest sum of distances to other points in the cluster:
  \vspace{-4pt}  $$\begin{array}{l} 
  D \{p\times p\} \rightarrow aclnnReduceSum \rightarrow sum\_dist\_temp \{k\times p\} \rightarrow \\aclnnArgMin \rightarrow min\_index\_temp \{1\}.
      \end{array}$$
    
    \item Updating the cluster's medoid index to the new one (executed on  Host):
  \vspace{-8pt}  $$\begin{array}{l}
  min\_ind\_temp \{1\} \rightarrow aclrtMemcpy \rightarrow min\_ind (C int);\\%[4pt]
    medoid[cl] = p\_ind[min\_ind],\mbox{where } cl \mbox{ is a current cluster.} 
      \vspace{-6pt}  \end{array}$$
      
    \item Steps 9--12 are performed for each cluster $cl = {1,\ldots,K}$.\\
    At the end of the cycle, the tensor $Xmed$ is updated:\\
$medoid (C array) \rightarrow aclrtMemcpy \rightarrow Xmed (aclTensor) \{k\times 1\}.$
    \item The conditions for exiting the loop are checked. If the $Xmed$ matrix is equal to the $P\_Xmed$ or $PP\_Xmed$ matrix, then this means that the optimal medoids for the clusters have been found and the $Xmed$ matrix no longer changes, the cycle ends. Otherwise, the steps 4--13 are repeated. Also, the algorithm terminates when %condition for exiting the loop is to reach 
    the maximum number of iterations has reached.
    \item The matrix $Xmed$ and point labels $L$ is uploaded to Host and the result is output.
 \vspace{-6pt}\end{enumerate}
 
  
\section{Computational experiments}
The computing environment used to execute and test our algorithm is equipped with a CPU Kunpeng-$920$ ($2.6$ GHz, aarch64), $1$ TB of RAM and $8$ NPU devices of Ascend $910$B with $64$GB of memory size. The operating system is EulerOS with the CANN version $8.0$ installed.

%\subsection{Datasets}
%In our experiments, two different types of datasets are used:
%\begin{itemize}
%    \item Syntetic randomly generated dataset
%    \item Second dataset
%\end{itemize}

The synthetic dataset is created by using a random data generator produ\-cing a uniform data distribution, so that no natural clustering structure exists. %This is to avoid biases introduced by pre-clustered data and represents the worst-case scenario for any clustering algorithm.


In Table \ref{tab:kmeans}, we report the results of running the NPU-accelerated $k$-means algorithm developed and its comparison against the two competing algorithms: classical $k$-means \cite{Lloyd82} running on CPU and Yinyang K-means \cite{DingZhaoShen15,TayGow21} implementing on GPU, which is one of the fastest $k$-means implementati\-ons. The column $m\_n\_k$ contains the number of data points in each instance ($m$), the number of data features ($n$),  the number of clusters to be found ($k$). The column $obj.val.$ shows the found objective value of the $k$-means problem. Next, for the three  algorithms, the number of iterations ($it$)  and the run time in seconds are specified.



%\subsection{Results for synthetic dataset}
%\vspace{-10pt}
{\small %\footnotesize
\begin{table}[htbp]
  \centering
  \caption{$k$-means algorithm}
    \vspace{-5pt} 
    \begin{tabular}{c|r|r|r|r|r|r|r}\hline
    \multirow{2}[2]{*}{\textbf{m\_n\_k}}  &\multirow{2}[2]{*}{\textbf{ obj.val $\:$}}  &\multicolumn{2}{c|}{\textbf{CPU}} &\multicolumn{2}{c|}{\textbf{GPU}}& \multicolumn{2}{c}{\textbf{NPU}}\\ 
    \cline{3-8}          
     &    & \multicolumn{1}{c|}{\textbf{it}} & \multicolumn{1}{c|}{\textbf{time, s}} & \multicolumn{1}{c|}{\textbf{it}} & \multicolumn{1}{c|}{\textbf{time, s}} & \multicolumn{1}{c|}{\textbf{it}} & \multicolumn{1}{c}{\textbf{time, s}}\\ \hline
   % \midrule
    100\_8\_10    & {3218.7}    & 8 & 0.0001 & 8     & 0.0012 & 8     & 0.7961\\ 
    1000\_16\_500   & {29366.9}   & 5 & 0.0579 & 5     & 0.0085 &5     & 0.7970\\
    10000\_64\_600   & {4060770.0} & 14 & 10.4625 & 14    & 0.2834 & 13    & 0.8509\\
    50000\_128\_800   & {4.73E+07}  & 33 & 363.6061 & 33    & 36.3916 & 26    & 1.3819 \\
    $\!\!$100000\_256\_1000$\!\!$  & {2.00E+08} & 55 & 3509.7585 & 55    & 292.8231 & 35    & 5.0792 \\
   $\!\!$500000\_1024\_1000$\!\!$  & 4.20E+09 & 65 & $\!\!$88672.9600 & \multicolumn{1}{r}{\footnotesize Out$\!\!$} & \multicolumn{1}{r|}{\footnotesize $\!\!\!$of memory}& 75 & $\!\!$29.1492\\ \hline
    \end{tabular}%
  \label{tab:kmeans}%
\end{table}%
 %\vspace{-6pt} 
}

On the small-scaled tests the NPU-accelerated $k$-means algorithm runs slower than both of its competitors. This behavior is not a surprise and comes from memory latency and data transfer overheads that makes algorithm slower. However, our NPU-accelerated implementation turns out to be $3000$ times faster than CPU $k$-means on the large-scale instances including $500\,000$ elements, $1024$ features and $1000$ clusters. In this test case, the Yinyang K-means on GPU fails due to the memory error. In the test with $m=100\,000$, the NPU-accelerated $k$-means is about $58$ times faster than the Yinyang \linebreak K-means on GPU which was announced to be the fastest implementation\linebreak of $k$-means.

Table \ref{tab:kmedoids} shows the results of testing the NPU-accelerated $k$-medoids algo\-rithm and the CPU-implemented conventional $k$-medoids. As well as for $k$-means, the NPU-accelerated $k$-medoids algorithm is slower on the small-scale tests up to 10000 items, 64 features, and 600 clusters. However, when clustering largest-scaled datasets, NPU-$k$-medoids algorithm turns out to be $445$ times faster than its CPU counterpart.

{\small 
 \begin{table}[htbp]
  \centering
  \caption{$k$-medoids algorithm}
  \vspace{-5pt} 
    \begin{tabular}{c|r|r|r|r|r}\hline
    \multirow{2}[2]{*}{\textbf{m\_n\_k}}  &\multirow{2}[2]{*}{\textbf{ obj.val $\:$}}  &\multicolumn{2}{c|}{\textbf{CPU}}& \multicolumn{2}{c}{\textbf{NPU}}\\ \cline{3-6}          
     &  &   \multicolumn{1}{c|}{\textbf{it}} & \multicolumn{1}{c|}{\textbf{time, s}} & \multicolumn{1}{c|}{\textbf{it}} & \multicolumn{1}{c}{\textbf{time, s}}\\ \hline
   % \midrule
    10000\_64\_600   & {6054170.0} & 2 & 1.3     & 2 & 3.33\\
    50000\_128\_800   & {7.42E+07}  & 2 & 22.2    & 2 & 3.41 \\
    100000\_256\_1000  & {3.32E+08}  & 2 & 124.8   & 2 & 4.11 \\
    500000\_1024\_1000  & {7.59E+09}  & 2 & 2882.38 & 2 & 6.47 \\ 
    \hline
    \end{tabular}%
  \label{tab:kmedoids}%
\end{table}%
}

 

  \vspace{-10pt}
\section{Conclusion}
In the paper, we proposed and developed NPU-accelerated implementati\-ons of the $k$-medoids and $k$-means algorithms. The computational study has demonstrated that the NPU-accelerated algorithms give a significant speedup compared to their sequential implementations as well as the known fastest alternative, Yinyang K-means \cite{DingZhaoShen15,TayGow21}.

%It is expected that due to their unique capabilities, NPUs are set to play a pivotal role in the future of AI, driving performance improvements and enabling new applications across various industries.

\bigskip
\begin{thebibliography}{10} % Ññûëêè íà ëèòåðàòóðó ïðèâîäÿòñÿ òîëüêî íà àíãëèéñêîì ÿçûêå (äàæå åñëè ó èñòî÷íèêà íåò àíãëîÿçû÷íîé âåðñèè).


\bibitem{VasUsh21} I. Vasilyev, A. Ushakov, \href{https://doi.org/10.1134/S1990478921040128}{\it  Discrete facility location in machine learning}, J. Appl. Ind.
Math.,  {\bf 15} (2021), 686--710.


\bibitem{VasUshGru15}
A.V. Ushakov, I.L. Vasilyev, T.V. Gruzdeva, \href{https://www.researchgate.net/publication/281698498_A_computational_comparison_of_the_p-median_clustering_and_k-means}{\it A computational comparison of the
p-median clustering and k-means}, Int. J. Artif. Intell., \textbf{13}:1 (2015), 229--242. 



\bibitem{MahNimVar12} M. Mahajan, P. Nimbhorkar, K. Varadarajan, \href{https://doi.org/10.1016/j.tcs.2010.05.034}{\it The planar k-means problem is NP-hard},
Theor. Comput. Sci., {\bf 442} (2012), 13--21.

\bibitem{AloDesHanPop09} D. Aloise, A. Deshpande, P. Hansen, P. Popat, \href{https://doi.org/10.1007/s10994-009-5103-0}{\it NP-hardness of Euclidean sum-of-squares
clustering}, Mach. Learn., {\bf  75} (2009), 245--248.

\bibitem{DolKel11} A. V. Dolgushev, A. V. Kelmanov, \href{https://doi.org/10.1134/S1990478911020050}{\it On the algorithmic complexity of a problem
in cluster analysis}, J. Appl. Ind. Math., \textbf{5}:2 (2011), 191--194.

\bibitem{KarivHakimi79}O. Kariv, S. Hakimi, \href{https://doi.org/10.1137/0137040}{\it An algorithmic approach to network location problems. I: The
p-centers}, SIAM J. Appl. Math., \textbf{37}:3 (1979), 513--538.

\bibitem{Megiddo79} N. Megiddo, \href{https://doi.org/10.1145/800133.804326}{\it Combinatorial optimization with rational objective functions}, Math.
Oper. Res., \textbf{4}:4 (1979), 414--424.

\bibitem{Lloyd82} S. Lloyd, \href{https://doi.org/10.1109/TIT.1982.1056489}{\it Least squares quantization in PCM}, IEEE Trans. Inf. Theory, \textbf{28}:2 (1982),
129--137.

\bibitem{SheLiuTsa17} X. Shen, W. Liu, I. W. Tsang, F. Shen, Q. S. Sun, 	\href{https://doi.org/10.1609/aaai.v31i1.10852}{\it Compressed k-means for large scale
clustering}, in: Proc. of the Thirty-First AAAI Conference on Artificial
Intelligence, February 4-9, 2017, San Francisco, California, USA., 2017, pp. 2527--
2533.


\bibitem{GongPawYan15} Y. Gong, M. Pawlowski, F. Yang, L. Brandy, L. Boundev, R. Fergus, \href{https://doi.org/10.1109/CVPR.2015.7298596}{\it Web scale photo
hash clustering on a single machine}, in: 2015 IEEE Conference on Computer Vision
and Pattern Recognition (CVPR), IEEE, New York, 2015, pp. 19--27.

\bibitem{DingZhaoShen15} Y. Ding, Y. Zhao, X. Shen, M. Musuvathi, T. Mytkowicz, \href{https://dl.acm.org/doi/10.5555/3045118.3045181}{\it Yinyang k-means: A drop-in
replacement of the classic k-means with consistent speedup}, in: Proc. of the 32nd International Conference on Machine Learning, {\bf 37}, ICML15, Lille, 2015, pp. 579--587.




\bibitem{YuChenChen20} Q. Yu, K.-H. Chen, J.-J. Chen,\href{https://doi.org/10.1007/978-3-030-60936-8\_23}{\it Using a set of triangle inequalities to accelerate
k-means clustering}, in: Similarity Search and Applications, LNCS, {\bf 12440}, Springer, Cham, 2020, pp. 297--311.





\bibitem{UshVas21} A. V. Ushakov, I. Vasilyev, \href{https://doi.org/10.1016/j.ins.2020.08.121}{\it Near-optimal large-scale k-medoids clustering}, Inf. Sci.,
{\bf 545} (2021), 344--362.

\bibitem{BayGhaHus18} M. Baydoun, H. Ghaziri, M. Al-Husseini, \href{https://doi.org/10.1007/s11227-018-2405-7}{\it CPU and GPU parallelized kernel k-means},
J. Supercomput., {\bf 74} (2018), 3975--3998.

\bibitem{BhiLeeMi15} J. Bhimani, M. Leeser, N. Mi, \href{https://doi.org/10.1109/HPEC.2015.7322467}{\it Accelerating k-means clustering with parallel
implementations and GPU computing}, in: 2015 IEEE High Performance Extreme
Computing Conference (HPEC), IEEE, New York, 2015, pp. 1--6.

\bibitem{CuoAngFarMar19} S. Cuomo, V. De Angelis, G. Farina, L. Marcellino, G. Toraldo,  \href{https://doi.org/10.1016/j.compeleceng.2017.12.002}{\it A GPU-accelerated
parallel k-means algorithm}, Comput. Elect. Eng., {\bf 75} (2019), 262--274.

\bibitem{DaoZouMez19} S. Daoudi, C. M. A. Zouaoui, M. C. El-Mezouar, N. Taleb,  \href{https://doi.org/10.1109/ICAEE47123.2019.9014783}{\it A comparative study of
parallel CPU/GPU implementations of the k-means algorithm}, in: 2019 International
Conference on Advanced Electrical Engineering (ICAEE), IEEE, New York, 2019,
pp. 1--5.

\bibitem{GaoWuLiaMen24} J. Gao, M.Wu, J. Liao, F. Meng, C. Chen,  \href{https://doi.org/10.1002/jcc.27470}{\it Clustering one million molecular structures
on GPU within seconds}, J. Comput. Chem., \textbf{45}:32 (2024), 2710--2718.

\bibitem{KruKra20} M. Krulivs, M. Kratochvil,  \href{https://doi.org/10.1145/3404397.3404426}{\it Detailed analysis and optimization of CUDA k-means
algorithm}, in: Proc. 49th International Conference on Parallel Processing, ACM, New
York, 2020, article 69.

%

\bibitem{LiFranPfa23} M. Li, E. Frank, B. Pfahringer, \href{https://doi.org/10.1007/s10618-022-00869-6}{\it Large scale k-means clustering using GPUs}, Data Min.
Knowl. Disc., {\bf 37} (2023), 67--109.

\bibitem{LiZhaoChuLiu13} Y. Li, K. Zhao, X. Chu, J. Liu, \href{https://doi.org/10.1016/j.jcss.2012.05.004}{\it Speeding up k-means algorithm by GPUs}, J. Comput.
Syst. Sci., \textbf{79}:2 (2013), 216--229.

\bibitem{LutBreRab18} C. Lutz, S. Breß, T. Rabl, S. Zeuch, V. Markl, \href{https://doi.org/10.1145/3211922.3211925}{\it Efficient k-means on GPUs}, in: Proc.
14th International Workshop on Data Management on New Hardware, DAMON'18,
ACM, New York, 2018, article 3.

\bibitem{MaLim13}
G. J. Lim, L. Ma,\href{https://doi.org/10.1016/j.cie.2012.10.008}{\it
GPU-based parallel vertex substitution algorithm for the p-median problem},
Computers \& Industrial Engineering, \textbf{64}:1
(2013),
381--388.

\bibitem{MoWangXhang22} Z. Mo, Y.Wang, Q. Zhang, G. Zhang, M. Guo, Y. Zhang, C. Shen, \href{https://doi.org/10.1007/978-3-031-15937-4_26}{\it The parallelization
and optimization of k-means algorithm based on MGPUSim}, in: 
Artificial Neural Networks and Machine Learning, ICANN 2022, Springer, Cham,
2022, pp. 309--320.

\bibitem{NelPal19} J. Nelson, R. Palmieri, \href{https://doi.org/10.1145/3303084.3309488}{\it Don't forget about synchronization! A case study of k-means on
GPU}, in: Proc. 10th InternationalWorkshop on Programming Models and Applications
for Multicores and Manycores, PMAM'19, ACM, New York, 2019, pp. 11--20.

\bibitem{TayGow21} C. Taylor, M. Gowanlock, \href{https://doi.org/10.1109/ICDE51399.2021.00163}{\it Accelerating the Yinyang k-means algorithm using the GPU},
in: 2021 IEEE 37th International Conference on Data Engineering (ICDE), IEEE,
New York, 2021, pp. 1835--1840.

%\bibitem{WuDingZhaiLiu24}
%S. Wu, Y. Ding, Y.  Zhai, J.   Liu,  J.  Huang, Z. Jian,  H. Dai, S. Di,   B. M. Wong, Z.  Chen,  F. Cappello, \href{https://doi.org/10.1109/CLUSTER59578.2024.00035}{\it FT K-Means: A High-Performance K-Means on GPU with Fault Tolerance}, in:
%2024 IEEE International Conference on Cluster Computing (CLUSTER), IEEE,
%New York,
%2024,
%pp. 322--334.

\bibitem{ParkJun09} H.-S. Park, C.-H. Jun, \href{https://doi.org/10.1016/j.eswa.2008.01.039}{\it A simple and fast algorithm for k-medoids clustering}, Expert.
Syst. Appl., \textbf{36}:2 (Part 2) (2009), 3336--3341.

\end{thebibliography}

%\bibliographystyle{splncs04}
%\bibliographystyle{elsart-num}
% \bibliography{bib-disser}

\end{document}
