# Detecting exact nearest neighbors with KMKNN

Aaron Lun1

1Cancer Research UK Cambridge Institute, Cambridge, United Kingdom

#### Revised: 27 September 2018

#### Package

BiocNeighbors 1.0.0

# 1 Introduction

The BiocNeighbors package provides an implementation of the k-means for k-nearest neighbors (KMKNN) algorithm, as described by Wang (2012). For a dataset with \(N\) points, the pre-training is done as follows:

1. Apply k-means clustering to all points, partitioning the data into \(\sqrt{N}\) clusters.
2. Compute the distance from each data point to its cluster center.
3. Store the cluster identities and distances.

For each query point, identification of the nearest neighbors is done as follows:

1. Start with a threshold distance \(d\) to the current kth-nearest neighbor (this can be set with arbitrary points).
2. Compute the distance from the query to each cluster center.
3. For any given cluster center, apply the triangle inequality on the query-center distance, the center-point distances and \(d\). Only compute query-point distances for points where the triangle inequality holds.
4. Update \(d\) with the new closest kth-nearest neighbor and repeat for the next cluster.

The pre-clustering arranges the points in a manner that effectively reduces the search space, even in high-dimensional data. Note that, while `kmeans` itself is random, the k-nearest neighbors result is fully deterministic11 Except in the presence of ties, see `?findKNN` for details..

The algorithm is implemented in a combination of R and C++, derived from code in cydar (Lun, Richard, and Marioni 2017). We observe 2-5-fold speed-ups in 20- to 50-dimensional data, compared to KD-trees in FNN and RANN (see https://github.com/LTLA/OkNN2018 for timings). This is consistent with results from Wang (2012).

# 2 Identifying k-nearest neighbors

The most obvious application is to perform a k-nearest neighbors search. We’ll mock up an example here with a hypercube of points, for which we want to identify the 10 nearest neighbors for each point.

    nobs <- 10000
    ndim <- 20
    data <- matrix(runif(nobs*ndim), ncol=ndim)

The `findKNN()` method expects a numeric matrix as input with data points as the rows and variables/dimensions as the columns. We indicate that we want to use the KMKNN algorithm by setting `BNPARAM=KmknnParam()` (which is also the default, so this is not strictly necessary here).

    fout <- findKNN(data, k=10, BNPARAM=KmknnParam())
    head(fout$index)

    ##      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
    ## [1,] 8797  137 9608 2387 1343 6746 3596 1590 8213  3373
    ## [2,] 8609 2956 3776  870 9502 8867 3195 7045 6876  7551
    ## [3,] 8891 5344 5020 5257 5073 1071 4841 7638  293  9643
    ## [4,] 8031 9606 8871 5307 9647 6906 1095   34 4624  6716
    ## [5,] 8637 6888 8570   55 9533 1184 2234 3769 5819  6400
    ## [6,] 8970 1312 4003 2625 8144 1746 5431 1130 3198   587

    head(fout$distance)

    ##           [,1]      [,2]      [,3]      [,4]      [,5]      [,6]      [,7]
    ## [1,] 0.8991621 0.9943289 0.9959295 1.0287310 1.0492102 1.0840995 1.0879738
    ## [2,] 0.8435266 0.9588740 0.9829271 1.0328775 1.0330807 1.0658921 1.0704479
    ## [3,] 0.8384262 0.9894006 0.9959119 1.0012585 1.0110796 1.0210744 1.0461125
    ## [4,] 0.7608781 0.8348439 0.8697829 0.8805191 0.8891034 0.9083789 0.9117551
    ## [5,] 0.6578826 0.9444255 0.9575260 0.9700629 0.9800220 1.0056250 1.0146012
    ## [6,] 0.9357125 0.9839811 1.0128718 1.0324356 1.0590379 1.0599838 1.0690042
    ##           [,8]      [,9]     [,10]
    ## [1,] 1.1050654 1.1130708 1.1136961
    ## [2,] 1.0907286 1.0918714 1.0949094
    ## [3,] 1.0674593 1.0716713 1.0859295
    ## [4,] 0.9298529 0.9368435 0.9485615
    ## [5,] 1.0333148 1.0415325 1.0484812
    ## [6,] 1.0734707 1.0824790 1.0833000

Each row of the `index` matrix corresponds to a point in `data` and contains the row indices in `data` that are its nearest neighbors. For example, the 3rd point in `data` has the following nearest neighbors:

    fout$index[3,]

    ##  [1] 8891 5344 5020 5257 5073 1071 4841 7638  293 9643

… with the following distances to those neighbors:

    fout$distance[3,]

    ##  [1] 0.8384262 0.9894006 0.9959119 1.0012585 1.0110796 1.0210744 1.0461125
    ##  [8] 1.0674593 1.0716713 1.0859295

Note that the reported neighbors are sorted by distance.

# 3 Querying k-nearest neighbors

Another application is to identify the k-nearest neighbors in one dataset based on query points in another dataset. Again, we mock up a small data set:

    nquery <- 1000
    ndim <- 20
    query <- matrix(runif(nquery*ndim), ncol=ndim)

We then use the `queryKNN()` function to identify the 5 nearest neighbors in `data` for each point in `query`.

    qout <- queryKNN(data, query, k=5, BNPARAM=KmknnParam())
    head(qout$index)

    ##      [,1] [,2] [,3] [,4] [,5]
    ## [1,] 6930 6113  245 6351 6875
    ## [2,] 6078 4456 5909  293 3715
    ## [3,] 6809 5827 2305 6650 9398
    ## [4,] 6314 6822 2176 2793 1932
    ## [5,] 5284 6789 7874 2918 1748
    ## [6,] 5210 2513  481  577 5914

    head(qout$distance)

    ##           [,1]      [,2]      [,3]      [,4]      [,5]
    ## [1,] 0.9449376 0.9863152 0.9921554 1.0143361 1.0228513
    ## [2,] 0.8024176 0.8813117 0.8853903 0.9374485 0.9379712
    ## [3,] 0.8659528 0.8883588 0.9239009 0.9443204 0.9604947
    ## [4,] 0.8654430 0.8874423 0.9023727 0.9682434 1.0236375
    ## [5,] 1.0825792 1.0833069 1.0837354 1.0885633 1.0926872
    ## [6,] 0.8119512 0.9247237 0.9389407 0.9727975 0.9748540

Each row of the `index` matrix contains the row indices in `data` that are the nearest neighbors of a point in `query`. For example, the 3rd point in `query` has the following nearest neighbors in `data`:

    qout$index[3,]

    ## [1] 6809 5827 2305 6650 9398

… with the following distances to those neighbors:

    qout$distance[3,]

    ## [1] 0.8659528 0.8883588 0.9239009 0.9443204 0.9604947

Again, the reported neighbors are sorted by distance.

# 4 Further options

Users can perform the search for a subset of query points using the `subset=` argument. This yields the same result as but is more efficient than performing the search for all points and subsetting the output.

    findKNN(data, k=5, subset=3:5)

    ## $index
    ##      [,1] [,2] [,3] [,4] [,5]
    ## [1,] 8891 5344 5020 5257 5073
    ## [2,] 8031 9606 8871 5307 9647
    ## [3,] 8637 6888 8570   55 9533
    ##
    ## $distance
    ##           [,1]      [,2]      [,3]      [,4]      [,5]
    ## [1,] 0.8384262 0.9894006 0.9959119 1.0012585 1.0110796
    ## [2,] 0.7608781 0.8348439 0.8697829 0.8805191 0.8891034
    ## [3,] 0.6578826 0.9444255 0.9575260 0.9700629 0.9800220

If only the indices are of interest, users can set `get.distance=FALSE` to avoid returning the matrix of distances. This will save some time and memory.

    names(findKNN(data, k=2, get.distance=FALSE))

    ## [1] "index"

It is also simple to speed up functions by parallelizing the calculations with the BiocParallel framework.

    out <- findKNN(data, k=10, BPPARAM=MulticoreParam(3))

For multiple queries to a constant `data`, the pre-clustering can be performed in a separate step with `buildNNIndex()`. The result can then be passed to multiple calls, avoiding the overhead of repeated clustering22 The algorithm type is automatically determined when `BNINDEX` is specified, so there is no need to also specify `BNPARAM` in the later functions..

    pre <- buildNNIndex(data, BNPARAM=KmknnParam())
    out1 <- findKNN(BNINDEX=pre, k=5)
    out2 <- queryKNN(BNINDEX=pre, query=query, k=2)

Advanced users may also be interested in the `raw.index=` argument, which returns indices directly to the precomputed object rather than to `data`. This may be useful inside package functions where it may be more convenient to work on a common precomputed object.

# 5 Session information

    sessionInfo()

    ## R version 3.5.1 Patched (2018-07-24 r75008)
    ## Platform: x86_64-w64-mingw32/x64 (64-bit)
    ## Running under: Windows Server 2012 R2 x64 (build 9600)
    ##
    ## Matrix products: default
    ##
    ## locale:
    ## [1] LC_COLLATE=C
    ## [2] LC_CTYPE=English_United States.1252
    ## [3] LC_MONETARY=English_United States.1252
    ## [4] LC_NUMERIC=C
    ## [5] LC_TIME=English_United States.1252
    ##
    ## attached base packages:
    ## [1] stats     graphics  grDevices utils     datasets  methods   base
    ##
    ## other attached packages:
    ## [1] BiocNeighbors_1.0.0 BiocParallel_1.16.0 knitr_1.20
    ## [4] BiocStyle_2.10.0
    ##
    ## loaded via a namespace (and not attached):
    ##  [1] Rcpp_0.12.19        bookdown_0.7        digest_0.6.18
    ##  [4] rprojroot_1.3-2     backports_1.1.2     stats4_3.5.1
    ##  [7] magrittr_1.5        evaluate_0.12       stringi_1.2.4
    ## [10] S4Vectors_0.20.0    rmarkdown_1.10      tools_3.5.1
    ## [13] stringr_1.3.1       parallel_3.5.1      xfun_0.4
    ## [16] yaml_2.2.0          compiler_3.5.1      BiocGenerics_0.28.0
    ## [19] BiocManager_1.30.3  htmltools_0.3.6

# References

Lun, A. T. L., A. C. Richard, and J. C. Marioni. 2017. “Testing for differential abundance in mass cytometry data.” Nat. Methods 14 (7):707–9.

Wang, X. 2012. “A Fast Exact k-Nearest Neighbors Algorithm for High Dimensional Search Using k-Means Clustering and Triangle Inequality.” Proc Int Jt Conf Neural Netw 43 (6):2351–8.