An Application of Storage-Optimal MatDot Codes for Coded Matrix Multiplication: Fast k-Nearest Neighbors Estimation

Utsav Sheth, Sanghamitra Dutta, Malhar Chaudhari, Haewon Jeong, Yaoqing Yang, Jukka Kohonen, Teemu Roos, Pulkit Grover

I INTRODUCTION

We consider the problem of finding the k nearest neighbors of a query point in a given high-dimensional dataset. To solve this problem efficiently, our goal is to speed up an existing algorithm by parallelizing it, and to make it resilient to stragglers . The k-nearest neighbor (k-NN) problem is often a first step used in a variety of real world applications including genomics , personalized search , network security , and web based recommendation systems .

In the era of Big Data, k-NN algorithms are often a bottleneck, as data and dimensionalities grow . There is a rich body of work on fast nearest neighbor retrieval. The existing work can be broadly classified into three categories. The first category of methods speeds up k-NN retrieval by reducing the space over which exact distance calculations are performed, by employing space partitioning data structures . The second category of techniques improves retrieval times by system level parallelism . Algorithms in the first and second categories are seldom scalable as they require that the whole dataset is held in (shared) memory for optimal performance. For large high-dimensional datasets, marshaling enough resources on a single system is challenging. The third category of methods aims to overcome this challenge by parallelizing storage and computation in a distributed setting. PANDA and DSI sharding are examples of data parallel implementations of distributed k-NN algorithms. While these techniques rely on specialized or high-performance hardware (e.g. Edison supercomputer for PANDA), the general trend in the distributed systems has been to use general purpose commodity systems . These approaches also do not address a more serious issue — the effects of node failures and slow nodes, or “stragglers”. Schroeder and Gibson observed as many as $1{,}159$ system failures per year at the Los Alamos National Laboratory. Dean et al. study a real Google service and observe that the slowest 5% of requests are responsible for half of the total 99th percentile latency.

Recently “Coded Computing” has been found to be very useful in combating stragglers and faults, by the efficient use of novel erasure-codes to create redundancy in computing. In this paper, we use one such coded computing technique called MatDot codes to speed up an approximate $k$ -NN algorithm called Multiple Random Projection Trees (MRPT) .

The problem of distributed matrix multiplication has always been of significant interest in the coded computing community . In our prior work , we proposed a novel coding technique called MatDot codes for the multiplication of two matrices (e.g. $\mathbf{X}^{T}\mathbf{Q}$ ) that outperform the recent Polynomial Codes with respect to recovery threshold, i.e., the number of worker nodes needed to wait for, out of the total nodes, under storage constraints. More specifically, for the problem of distributed matrix multiplication under the storage constraint that each node can only store a fixed $\frac{1}{m}$ fraction of each matrix, MatDot codes achieve the optimal recovery threshold of $2m-1$ as compared to Polynomial codes that have a recovery threshold of $m^{2}$ .

We note that for matrix multiplication $\mathbf{X}^{T}\mathbf{Q}$ or matrix-vector product $\mathbf{X}^{T}\mathbf{q}$ under the constraint that only $\frac{1}{m}$ fraction of each operand can be stored at each node, MatDot codes use vertical block-partitioning of the first matrix $\mathbf{X}^{T}$ as compared to other existing strategies that use either horizontal partition or a combination of both (see for discussion on partitioning). Interestingly, as it turns out, the vertical partitioning of the first matrix is better suited for the coded MRPT problem formulation. This is because we are required to perform the multiplication $\mathbf{X}(S)^{T}\mathbf{Q}$ or $\mathbf{X}(S)^{T}\mathbf{q}$ in a distributed fashion, where $\mathbf{X}(S)^{T}$ denotes a sub-matrix of $\mathbf{X}^{T}$ consisting of the rows of $\mathbf{X}^{T}$ indexed in a set $S$ , such that, $\bm{X}^{T}$ is known in advance but the set $S$ becomes available only in the online phase (real-time).

This coded computing problem ensues from our broader goal in this work, which is to speed up MRPT by employing data and model parallelism and straggler-tolerant computing techniques. We differentiate here between data parallelism and model parallelism. In data parallelism, different nodes process different pieces of data, but each node performs all the computations relevant to the entire model. In model parallelism, the model itself is parallelized across nodes.

MRPT partitions the search space to retrieve approximate k nearest neighbors of the query $\mathbf{q}$ . It uses a combination of random projection trees and voting to achieve fast queries and high accuracy. In this paper, we propose two enhancements to the MRPT algorithm by parallelizing it in a distributed setting. Our contributions are as follows:

We propose a distributed implementation of MRPT exploiting data parallelism that experimentally demonstrates faster queries than a single node implementation, even when using CPUs with lower clock speeds. Additionally, the cloud-based virtual machines we use in our experiments for parallel MRPT have a non-zero steal time, i.e., they may be required to wait while others are being served.

We formulate a coded computing problem for MRPT, and then apply coded matrix multiplication strategies, namely MatDot and Systematic MatDot codes to further reduce the query time for the model parallel architecture in a system that is prone to straggling.

The rest of the paper is organized as follows. In Section II we explain the MRPT algorithm and describe how it reduces the search-space through projections and voting. In Section III we introduce the Data Parallel Model Implementation of MRPT. In Section IV we introduce our proposed model parallel implementation of MRPT and then describe the application of MatDot codes and systematic MatDot codes in our model parallel architecture to achieve a lower recovery threshold under straggling. The model parallel architecture is ideal for applications where system components are unreliable and accuracy is important. In Section V we experimentally demonstrate the advantages of our approach. A conclusion is provided in Section VI.

II Preliminaries

This section briefly describes the two stages of the MRPT algorithm: (i) off-line index construction stage and (ii) on-line query stage. Assume that we are given a $d$ -dimensional dataset $\mathcal{X}$ consisting of $N$ points, represented as a $d\times N$ matrix $\mathbf{X}$ . Given a query point $\mathbf{q}$ , the problem of k-nearest neighbors involves finding a set of points $\kappa\subseteq\mathcal{X}$ such that $|\kappa|=k$ and $dist(\mathbf{x},\mathbf{q})\leq dist(\mathbf{y},\mathbf{q})$ for each $\mathbf{x}\in\kappa$ , $\mathbf{y}\in\mathcal{X}\backslash\kappa$ , and the function $dist(\cdot)$ is the distance function in the d-dimensional Euclidean space given by:

where $\mathbf{u}$ and $\mathbf{v}$ are two vectors in this space.

In the MRPT algorithm, a sparse $d$ -dimensional random projection vector $\mathbf{r}$ is chosen, in which each entry $r_{i}$ is sampled from the following distribution:

Given a $d$ -dimensional query vector $\mathbf{q}$ , the first step in the MRPT query stage is to generate a candidate set of indices (pruned data-point indices) $S\subset\{1,2,\ldots,N\}$ such that $|S|$ $\ll N$ .

For each Random Projection (RP) tree $t\in$ $\mathcal{T}$ , at each level the query vector $\mathbf{q}$ is projected onto the random vector $\mathbf{r}$ for that level and then assigned a branch based on whether its value is greater than or less than the median of the projections of all other data-points with $\mathbf{r}$ . This process is then repeated recursively until a leaf is reached.

Each tree had already partitioned the dataset $\mathcal{X}$ into $2^{l}$ cells or leaves. For $1\leq t\leq T$ , let $f_{t}(\cdot)$ be defined as:

The candidate set of indices (pruned points) $S$ can then be finally chosen as follows:

Here, $\nu$ is a pre-configured parameter known as the voting threshold. Thus, the set $S$ denotes the set of indices $\subset\{1,2,\ldots,N\}$ for which at least $\nu$ trees have found the corresponding data-point $\mathbf{x}_{j}$ in the same cell as $\mathbf{q}$ .

II-B Coded Matrix Multiplication

Coded computing combines distributed numerical algorithms and error correcting codes (ECCs) to mitigate unreliable processors and randomness in their response time. In this work, we focus on coded matrix multiplication as matrix multiplication is the main bottleneck in MRPT algorithm. Coded computing has been used extensively for different computation objectives such as neural-network training, FFT, iterative computing, distributed regression, and convolutions (see for review).

System Model: We want to compute $\bm{C}=\bm{A}\bm{B}$ where $\bm{A}$ and $\bm{B}$ are $N$ -by- $N$ matrices. A master node distributes the computation to $P$ worker nodes. A worker node has limited memory/computing power, so each node can receive the $1/m$ -th fraction of matrices $\bm{A}$ and $\bm{B}$ . After completing its computation, a worker reports the result to a fusion node. Recovery threshold is defined as the worst-case number of workers needed to recover the final result.

Let $K$ be the number of workers needed to complete the computation if all worker nodes are reliable ( $K$ is different depending on how we split the matrix). In reality, some processors are significantly slower than the others due to queuing delays or random faults in the processor. Without any reliability measure to alleviate straggler problems, computation completion time would be dominated by few stragglers. Our aim is to use more than $K$ worker nodes by adding some redundancies, which are carefully designed by applying the ideas from coding theory, so that the whole computation can be resilient to stragglers.

MatDot codes: The matrix $\mathbf{A}$ is split vertically into $m$ column blocks, and $\mathbf{B}$ is split horizontally into $m$ row blocks:

where $\mathbf{A}_{i},\mathbf{B}_{i}$ ( $i=1,\ldots,m$ ) are $N\times N/m$ and $N/m\times N$ dimensional submatrices, respectively.

The matrices $\mathbf{A}$ and $\mathbf{B}$ are then encoded as polynomials:

A master node distributes encoded matrices, $p_{\mathbf{A}}(\alpha_{i})$ and $p_{\mathbf{B}}(\alpha_{i})$ to the $i$ -th worker node ( $i=1,\ldots,P$ ). Then the $i$ -th worker node computes the following product at $x=\alpha_{i}$ :

and returns the result to the master node. Note that the coefficient of $x^{m-1}$ in $p_{\mathbf{C}}(x)$ is $\mathbf{C}=\sum_{i=1}^{m}A_{i}B_{i}$ . Since $p_{\mathbf{C}}(x)$ is a polynomial of degree $2m-2$ , its coefficients can be recovered by the master node as soon as it receives the values of $p_{\mathbf{C}}(x)$ at any $2m-1$ distinct points. Hence the recovery threshold is $K=2m-1$ . This is provably the optimal recovery threshold, when a worker node can store $\frac{1}{m}$ -th fraction of each input matrix .

Systematic MatDot codes: A code is called systematic if, for the first $m$ worker nodes, the output of the $r$ -th worker node is the product $\mathbf{A}_{r}\mathbf{B}_{r}$ . We refer to the first $m$ worker nodes as systematic worker nodes. Having systematic nodes is useful because if all the systematic nodes complete their computation in time, there is no need for decoding. Systematic MatDot codes are achieved by applying different encoding polynomials. Let $p_{\mathbf{A}}(x)=\sum_{i=1}^{m}\mathbf{A}_{i}L_{i}(x)$ and $p_{\mathbf{B}}(x)=\sum_{i=1}^{m}\mathbf{B}_{i}L_{i}(x)$ where $L_{i}(x)$ is defined as follows for $i\in\{1,\ldots,m\}$ :

Using these polynomials, the worst-case recovery threshold remains the same as non-systematic MatDot codes .

MatDot and Systematic MatDot codes use vertical partitioning of the first matrix which is well suited for the problem of coded MRPT as compared to strategies that use horizontal partitioning or a combination of horizontal and vertical partitioning for the first matrix. This is because in the coded MRPT formulation, we are required to perform the multiplication $\mathbf{X}(S)^{T}\mathbf{Q}$ or $\mathbf{X}(S)^{T}\mathbf{q}$ in a distributed fashion, where $\mathbf{X}(S)^{T}$ denotes a sub-matrix of $\mathbf{X}^{T}$ consisting of the rows of $\mathbf{X}^{T}$ indexed in a set $S$ , such that the set $S$ is available only in the online phase (real-time).

III DATA PARALLEL MRPT

Now, we model MRPT as a problem in data parallelism and describe our first strategy to parallelize the algorithm. Consider a distributed computing cluster having a single master node and $P$ worker nodes as shown in Fig. 1.

Given a $d\times N$ matrix X representing the set of data-points $\mathcal{X}$ and a cluster with $P$ worker nodes, we randomly split $\mathbf{X}$ vertically into $P$ disjoint, vertical partitions $\mathbf{X}_{i}$ for $i\in{\{1,2,\ldots,P\}}$ .Thus, $\mathbf{X}=\begin{bmatrix}\mathbf{X}_{1}|\mathbf{X}_{2}|\dots|\mathbf{X}_{P}\end{bmatrix},$ where each $\mathbf{X}_{i}$ is of dimension $d\times\frac{N}{P}.$

We distribute each partition $\mathbf{X}_{i}$ across the $P$ worker nodes in the cluster such that the worker $W_{i}$ contains partition $\mathbf{X}_{i}$ . Each worker then runs the MRPT index construction algorithm described in Section II-A to build a local MRPT set of trees $\mathcal{T}_{i}$ from its partition $\mathbf{X}_{i}$ of $\mathbf{X}$ .

Given a query $\mathbf{q}$ , the master node transmits $\mathbf{q}$ to each worker node $W_{i}$ . Then $W_{i}$ uses the MRPT query algorithm described in Section II-A on its trees $\mathcal{T}_{i}$ to determine the $\tau$ nearest neighbors of $\mathbf{q}$ . We use the same voting threshold $\nu$ in all workers. After local voting, a local exact distance calculation step is performed to narrow down to $\tau$ data-points. Finally, each worker node $W_{i}$ then returns the indices of its set of $\tau\;(\geq k)$ nearest neighbors $N_{i}$ , from the partition $\mathbf{X}_{i}$ of $\mathbf{X}$ , to the master node along with their exact distances from the query $\mathbf{q}$ . Then, the master node determines the final set of k nearest neighbors to $\mathbf{q}$ from all the $P\tau$ candidate data-point indices received from all the worker nodes, i.e., the indices of all the data-points in the set $\cup_{i=1}^{P}N_{i}$ , by sorting the data-points based on their exact distances.

: Recall that $\kappa\subseteq\mathcal{X}$ (with $|\kappa|=k$ ) is the set of the true k nearest neighbors of $\mathbf{q}$ . We let $\kappa_{i}\subseteq\mathbf{\kappa}$ be the set of the true nearest neighbors to $\mathbf{q}$ that lie in worker $W_{i}$ . To achieve maximum accuracy, we must have $\kappa_{i}\subseteq N_{i}$ and therefore the value of $\nu$ chosen for the system must be such that $|<math xmlns="http://www.w3.org/1998/Math/MathML" display="block"><semantics><mrow><msub><mi>κ</mi><mi>i</mi></msub></mrow><annotation encoding="application/x-tex">\kappa_{i}</annotation></semantics></math>κi|$ is much less than $\tau$ . In fact, a higher value of $\tau$ implies higher chance of containing the desired set $\kappa_{i}$ , though it comes with increased communication cost from worker to master and increased computation at the master node.

IV MODEL PARALLEL MRPT

In this section we discuss a model parallel architecture for approximate k nearest neighbor search using MRPT. We then propose two enhancements to the model parallel architecture that apply coded distributed matrix multiplication techniques that achieve the optimal recovery threshold in a system that is prone to straggling.

Consider the d-dimensional dataset $\mathcal{X}$ as before. In the model parallel architecture, given a query $\mathbf{q}$ , we first find the possible candidate set of indices (pruned indices) S using the recursive algorithm described in Section II-A. Now the search space for the true nearest neighbors $\kappa$ reduces to the set of data-points whose indices are in $S$ , i.e., $\kappa\subseteq\{\mathbf{x}_{j}:j\in S\}$ .

To find the set $\kappa$ , we compute the exact Euclidean distance from each data-point $\mathbf{x}_{j}$ (for $j\in S$ ) to the query point $\mathbf{q}$ . Examining the terms constituting the Euclidean distance in (1), the Euclidean norm $\|<math xmlns="http://www.w3.org/1998/Math/MathML" display="block"><semantics><mrow><msub><mi mathvariant="bold">x</mi><mi>j</mi></msub></mrow><annotation encoding="application/x-tex">\mathbf{x}_{j}</annotation></semantics></math>xj\|$ for each $\mathbf{x}_{j}\in\mathcal{X}$ can be precomputed and the same can be done to get $\|<math xmlns="http://www.w3.org/1998/Math/MathML" display="block"><semantics><mrow><mi mathvariant="bold">q</mi></mrow><annotation encoding="application/x-tex">\mathbf{q}</annotation></semantics></math>q\|$ . We must now only compute the dot product $\mathbf{x}_{j}$ $\cdot$ $\mathbf{q}$ to obtain the Euclidean distances from each $\mathbf{x}_{j}$ to $\mathbf{q}$ .

To do this, we first represent the data-points indexed in the set S as a $d\times|S|$ matrix $\mathbf{X}(S)$ that contains only the data-points $\mathbf{x}_{j}$ (columns of $\mathbf{X}$ ) such that $j\in S$ . The transpose of this matrix is the $|S|\times d$ matrix $\mathbf{X}(S)^{T}$ that essentially denotes all the rows of the matrix $\mathbf{X}^{T}$ indexed in $S$ .

Consider the column vector $\mathbf{w}$ such that $\mathbf{w}=\mathbf{X}(S)^{T}\mathbf{q}$ . Note that each element of the vector $\mathbf{w}$ corresponds to the dot-product $\mathbf{x}_{j}\cdot\mathbf{q}=\mathbf{x}_{j}^{T}\mathbf{q}$ for some $j\in S$ . The Euclidean distance from $\mathbf{x}_{j}$ to $\mathbf{q}$ can now be determined as all the terms in (1) are known to us, which includes the individual norms as well as the dot product $\mathbf{x}_{j}\cdot\mathbf{q}$ .

The problem thus reduces to the following: compute the vector $\mathbf{w}=\mathbf{X}(S)^{T}\mathbf{q}$ in a distributed computing cluster where $\mathbf{X}^{T}$ is known in advance but the set of indices $S$ become available only in the online phase (real-time). Since the computation $\mathbf{w}=\mathbf{X}(S)^{T}\mathbf{q}$ is the only stage of the algorithm that must be done at runtime (in the online stage) and scales linearly with $d$ , we now discuss several scalable strategies that compute vector $\mathbf{w}$ in a distributed setting.

IV-B Uncoded Distributed Matrix-Vector Multiplication

We split $\mathbf{X}^{T}$ into $P$ equal partitions as follows (see Fig. 2):

Now consider a cluster consisting of one master node and $P$ worker nodes as shown in Fig. 3. Each partition $\mathbf{X}^{T}_{i}$ is distributed across the worker nodes such that worker $\textit{W}_{i}$ stores the partition ${\mathbf{X}_{i}}^{T}$ in advance (off-line).

Note that, if $\mathbf{X}^{T}$ is partitioned using the strategy just discussed, the matrix $\mathbf{X}(S)^{T}$ also gets partitioned as follows:

In the online phase, we only split the query $\mathbf{q}$ into $P$ equal partitions, $\{\mathbf{q}_{i}\,:\,i\in\{1,\ldots,P\}\}$ (again see Fig. 2). The product $\mathbf{w}$ can then be expressed as:

Given a query $\mathbf{q}$ for which the k nearest neighbors must be determined, the master node first computes the possible candidate set S for $\mathbf{q}$ from its MRPT index set of trees $\mathcal{T}$ and then transmits the set S and partition $\mathbf{q}_{i}$ of $\mathbf{q}$ to worker node $\textit{W}_{i}$ . For every S, each worker node $\textit{W}_{i}$ only fetches the matrix $\mathbf{X}_{i}(S)^{T}$ from $\mathbf{X}^{T}_{i}$ already stored in its memory. It then computes the product $\mathbf{X}_{i}(S)^{T}\mathbf{q}_{i}$ and returns the resulting vector to the master node. The master node can thus compute the vector $\mathbf{w}$ by adding the results using (11), and determine the $k$ nearest neighbors using the exact distances.

IV-C Coded Distributed Matrix-Vector Multiplication using MatDot Codes

In order to successfully compute the vector $\mathbf{w}$ using the strategy described in Section IV-B, the master node must wait for every worker node $\textit{W}_{i}$ to successfully return the product $\mathbf{X}_{i}(S)^{T}\mathbf{q}_{i}$ . In a straggler-prone environment, this might cause unprecedented delays in computation. Thus, to avoid waiting for all nodes and be able to recover the matrix-vector product by only waiting for some out of all workers to finish, we will now apply the MatDot-based distributed matrix multiplication strategy .

We partition the matrix $\mathbf{X}^{T}$ vertically again, but into $m$ partitions instead of $P$ as follows:

We then use the following encoding polynomial:

The rows of $P_{\mathbf{X}^{T}}(\beta)$ indexed in set $S$ actually represent the following polynomial:

We will be referring to this observation later.

Now, given a cluster with a master node and $P$ worker nodes, as shown in Fig. 4, each worker node $\textit{W}_{i}$ is initialized with a different $\beta_{i}$ , using which it computes the polynomial $\textit{P}_{\mathbf{X}^{T}}$ ( $\beta_{i}$ ) in (13). This encoding step can be performed off-line as $\mathbf{X}^{T}$ is known in advance.

During the online stage, given a query $\mathbf{q}$ for which the k nearest neighbors must be determined, the master node first partitions $\mathbf{q}$ into $m$ parts: $\{\mathbf{q}_{j}:\textit{j}\in\{1,2,\ldots,\textit{m}\}\}$ . We then use the following encoding polynomial:

As in Section IV-B, the master node first determines the candidate set S. It then transmits S and the encoded query $P_{\mathbf{q}}(\beta_{i})$ obtained from (15) to worker $\textit{W}_{i}$ . The worker $\textit{W}_{i}$ then fetches only the matrix $P_{\mathbf{X}(S)^{T}}(\beta_{i})$ from its stored $P_{\mathbf{X}^{T}}(\beta_{i})$ (recall (14)) which essentially denotes all the rows of $P_{\mathbf{X}^{T}}(\beta_{i})$ indexed in $S$ . Then, it computes the product $P_{\mathbf{X}(S)^{T}}(\beta_{i})P_{\mathbf{q}}(\beta_{i})$ and returns the result to the master node.

The coefficient of $\beta^{m-1}$ in the polynomial $P_{\mathbf{X}(S)^{T}}(\beta)P_{\mathbf{q}}(\beta)$ turns out to be our desired desired matrix-vector product $\mathbf{X}(S)^{T}\mathbf{q}=\sum_{j=1}^{m}\mathbf{X}_{j}(S)^{T}\mathbf{q}_{j}$ from the property of MatDot codes. We need to evaluate the polynomial at only $2m-1$ distinct points so as to determine the coefficient for every power of $\beta$ . The master node must therefore wait for at least $2m-1$ worker nodes following which it can determine the term $\sum_{i=1}^{m}\mathbf{X}_{i}(S)^{T}\mathbf{q}_{i}$ using polynomial interpolation. We then follow the strategy of comparing the exact distances in Section IV-A to obtain the set of the k nearest neighbors to $\mathbf{q}$ .

IV-D Coded Distributed Matrix-Vector Multiplication using Systematic MatDot Codes

In our prior work where we proposed MatDot Codes, we also introduced their systematic variant. Their advantage is that, while for MatDot Codes the recovery threshold is always $2m-1$ , for systematic MatDot Codes one might sometimes only need $m$ nodes to finish, although $2m-1$ is the worst-case value. In this section, we apply the systematic MatDot code to the MRPT problem.

Similar to the previous case, we first partition $\mathbf{X}^{T}$ vertically into $m$ partitions, but then use a different encoding function:

Consider a cluster consisting of one master node and $P$ worker nodes. Worker $\textit{W}_{i}$ is assigned a value $\textit{$ \beta $}_{i}$ using which it computes the polynomial in (16). This encoding is performed in advance, in the off-line stage. Interestingly, the workers $\{W_{i}:i=1,2,\ldots,m\}$ turn out to be the systematic worker nodes, which contain uncoded partitions of $\mathbf{X}^{T}$ .

In the online phase, given a query $\mathbf{q}$ , the master node first partitions $\mathbf{q}$ into $m$ partitions and then uses the following encoding function:

where $L_{j}(\beta)$ is given by (17). The master node then transmits candidate set $S$ and $P_{\mathbf{q}}(\beta_{i})$ to each worker node. Worker $\textit{W}_{i}$ is responsible for computing the product $P_{\mathbf{X}(S)^{T}}(\beta_{i})P_{\mathbf{q}}(\beta_{i})$ . We first consider the case when the first m workers to successfully complete their computation are the systematic worker nodes. We can then obtain the vector $\mathbf{w}$ as follows:

If the results of the first $2m-1$ successful workers do not contain results from the $m$ systematic nodes, then the master interpolates the polynomial $P_{\mathbf{X}(S)^{T}}(\beta)P_{\mathbf{q}}(\beta)$ . It then computes this polynomial product at each $\beta_{i}$ $\in$ { $\beta_{1}$ , $\beta_{2}$ ,…, $\beta_{m}$ }. Finally, it computes the vector $\mathbf{w}$ using (19). Note that in the ideal case, i.e., when all the systematic worker nodes finish first, we only needed $m$ nodes to finish as opposed to the worst-case recovery threshold of $2m-1$ . We can now proceed with the steps of comparing exact distances (see Section IV-A) to retrieve the k nearest neighbors of $\mathbf{q}$ .

V EXPERIMENTAL RESULTS

In this section, we evaluate the effectiveness of data and model parallel MRPT in terms of both accuracy and speed. All of our experiments were conducted on Amazon Elastic Compute Cloud instances .

The STL-10 dataset is a dataset of $N=100000$ images each of dimension $d=9216$ used in unsupervised image classification algorithms, while GIST is a popular dataset with $N=1000000$ and $d=960$ used in ANN algorithms. These datasets provide us with a good mix of dimensionality and number of datapoints to evaluate our proposed strategies. The MRPT parameters used for the experiments are provided in Table I.

We evaluate the accuracy of our implementation using recall defined as: $\dfrac{|\kappa_{MRPT}\cap\kappa|}{k}$ , where $\mathbf{q}$ is the query whose true k nearest neighbors is the set $\kappa$ and ${\kappa_{MRPT}}$ is the set of k nearest neighbors returned by the algorithm.

For the single node MRPT baseline, we used a compute optimized c5.large instance with two 3GHz Intel Xeon Platinum processors and 4 GB of memory. For the data parallel and model parallel architectures, we used t2.medium instances with two 2.3 GHz Intel Broadwell processors and 4 GB of memory. All instances are provisioned with 40GB HDD secondary storage. Note that the hardware used for our baseline experiments is superior to that used to evaluate our parallel strategies. In all experiments we ran 500 queries sequentially with $k=10$ . We consider the average result of 50 runs for each experiment. The experiments were conducted in a cluster consisting of 1 master node and 16 worker nodes.

In the experiments with MatDot codes and systematic MatDot codes, we used encoding polynomials of degree 2.

All our experiments are conducted on systems with limited memory. Thus the MRPT algorithm has to use the disk to hold the index $\mathcal{T}$ and the actual data points. In comparison, the data parallel and model parallel architectures use less memory and avoid disk penalties. In the data parallel architecture, we reduce the amount of points each worker must hold by a factor of P. Additionally, the MRPT indexes become smaller as each worker holds a smaller fraction of the dataset. In the model parallel architecture, the master node can discard the components for each point after tree construction; it only needs the index. At each worker node, the memory requirement is at least halved.

V-B Simulating Stragglers on Amazon Web Services

To demonstrate the effects of stragglers in the model parallel architecture, we sample the minimum time $T_{i}$ a worker $W_{i}$ must take to complete a matrix multiplication for a query from the shifted Exponential and Weibull distributions shown in (20) and (21) respectively.

Here, $l_{i}$ is the number of row vectors loaded at worker $W_{i}$ for matrix multiplication, $a_{i}>0$ is the shift parameter, $\alpha_{i}>0$ is the shape parameter for the Weibull distribution, $\mu_{i}>0$ is the straggling parameter for $W_{i}$ , and $t\geq{a_{i}}{l_{i}}$ . For our experiments, we set each $\mu_{i}$ , $a_{i}$ , and $\alpha_{i}$ to some constants $\mu$ , $a$ , and $\alpha$ to maintain homogeneity in minimum computation times across workers. The parameters chosen are as follows: Shifted Exponential ( $a=0.0000001$ , $\mu=15$ ) and Weibull ( $a=0.2$ , $\mu=2$ , $\alpha=0.5$ ) for both the datasets. Note that, for model parallel MRPT, $l_{i}=|S|$ . We use a similar strategy to simulate straggling in data parallel MRPT with $l_{i}$ set to the average of $|S|$ for the set of test queries.

V-C Results

Our experimental results are provided in Tables II and also illustrated in Fig. 5. For completion, we also include our obtained recall values here: STL-10: Data Parallel (0.9632), Model Parallel (0.9648). GIST: Data Parallel (0.9430), Model Parallel (0.9350).

In all our experiments, both data and model parallel MRPT outperform the single node implementation. As shown in Fig. 5(a) and Fig. 5(b), data parallel MRPT is significantly faster than single node MRPT. This is due to the smaller size of the MRPT index and absence of disk penalties at workers. It can be seen that the data parallel strategy has better performance when compared to uncoded model parallel strategy because of its embarrassingly parallel design. Owing to this design, the data parallel strategy could scale linearly with the number of nodes. However, these scaling benefits in query execution come at a cost, as the random projection trees computed at each node do not contain all the data points and hence the candidate set generated by each node could contain lesser true positives. To offset this condition, we might have to lower the voting threshold in-order to generate a better candidate, while causing more communication overheads and hence lower query execution time. The note mentioned in Section III explains this case in more detail. For our experiments, we do not lower the voting threshold as the loss in recall for STL-10 and GIST is not significant.

The model parallel architecture results in a high communication cost as the candidate set has to be transmitted to each worker node. However, we find that if the MRPT parameters for a dataset are sufficiently tuned, Algorithm 1 will generate a smaller candidate set over which exact distance calculations must be performed. For our experiments, the algorithm was able to reduce the search space for STL-10 from 100000 datapoints to 3025 and for GIST from 1000000 to 14934 datapoints on average per query. Additionally, the model parallel strategy does not suffer from the accuracy related issues as compared to the data parallel architecture as the candidate set is generated using the entire dataset and not parts of it separately.

Both the strategies outperform baseline single node MRPT despite running on inferior hardware. The parallel strategies are not limited by memory constraints as in the case of single node MRPT. These strategies are therefore very useful when large datasets do not fit in memory. The coded model parallel strategy also makes the algorithm tolerant to slow nodes and failures in a distributed setting.

Fig. 5(c) and Fig. 5(d) show that model parallel MRPT outperforms single node MRPT. They also outperform data parallel MRPT under simulated straggling. This is of significance to real world systems where straggling may manifest as unreliable nodes or network delays. Fig. 5(e) and Fig. 5(f) show the benefits of coded matrix multiplication as opposed to the uncoded model parallel architecture under simulated straggling. Both MatDot codes and systematic MatDot codes are consistently faster than the uncoded approach. Fig. 5(e) also shows that systematic MatDot codes is able to outperform MatDot codes owing to its lower recovery threshold.

VI CONCLUSIONS

We proposed two approaches to parallelize the MRPT algorithm in a distributed setting. We also applied the MatDot code based distributed matrix multiplication strategy to reduce the recovery threshold in a system that is prone to stragglers. We showed that our parallelization strategies can achieve faster queries than the single node MRPT algorithm under limited memory. Our results demonstrate the benefits of applying MatDot code and systematic MatDot code to the model parallel architecture in a system with simulated stragglers. In our experiments we observed large floating point errors when inverting high degree Vandermonde matrices for polynomial interpolation. As future work, we will experiment with strategies to reduce the condition number of a Vandermonde matrix so that we can employ polynomials of higher degrees, and thus apply the MatDot code and systematic MatDot code with a larger number of worker nodes. Another possibility is to perform the computations in exact rational arithmetic. This would eliminate rounding errors, but the effect on runtime needs to be analyzed.

Acknowledgements

This work was supported by NSF CNS-1702694, the Academy of Finland under the WiFIUS program, the Academy of Finland COIN CoE and NSF CCF 1350314.