Fastest array-by-array nearest neighbour search Intro -
Nearest neighbour comparison between two arrays is something that can be done very efficiently when only 1 variable per array is involved, because it can be done by sorting and then performing a rolling comparison. However, for a 2 variable comparison this sort operation cannot be done anymore and another smart way has to be thought of. Hence the challenge:

Challenge -
The challenge is to write a code that can, for each row in array A, find a row in array B such that (A[i,1]-B[j,1])^2+(A[i,2]-B[j,2])^2 is minimal and that then stores the corresponding row number j. The perfect minimum (0) is a perfect match of A[i,1] with B[j,1] and of A[i,2] with B[j,2].

As a code example I have a very clunky loop written in R

N=nrow(SetA)
for (i in 1:N){

      #calculate vector with distance of row i to the entire SetB
      vec=with(SetB,(SetA[i,x]-x)^2+(SetA[i,y]-y)^2)

      #calculate minimum of that vector and store the index
      indexVec[i,ind:=which.min(vec)]
}


Some sidenotes -
If it matters (not sure it will) assume that both arrays have 10k rows. Ties can be decided at 'random' e.g. determined by processing sequence or whatever way is more convenient in your code.

Example in/out -
For clarity, here is the sample input

SetA=[1 3;2 4;3 5;4 6;5 7]
SetB=[2 2.5; 4 3.1; 2 2.0; 3 3.0; 6 5.0;]


and the expected sample output

indexVec=[1;4;4;5;5]


where the values in indexVec indicate where the nearest neighbour of row i in SetA is located in SetB.

Fastest code (as tested on my i5-4210U) wins. The testset will consist of 2 arrays of 10k rows and 2 columns filled with floating point numbers between -1E10 and +1E10. Of course seeded with the same seed to remove any unfairness.

Feel free to use whatever language, but I am also very interested in an optimized version of the R code above.
C++, .0007 seconds

Uses a simple grid to place all the points of B, then looks up the points of A starting at the grid point to which it maps.

#include &lt;stdio.h&gt;
#include &lt;stdlib.h&gt;
#include &lt;thread&gt;
#include &lt;sys/time.h&gt;

double mytime() {
  struct timeval tv;
  gettimeofday(&amp;tv, NULL);
  return tv.tv_sec + tv.tv_usec * 1e-6;
}

typedef float vector[2];

#define T 4
// N must be a multiple of T
#define N 10000

// Width of grid.  Should be about sqrt(N)
#define W 100

vector a[N];
vector b[N];
int min[N];

float dist(vector x, vector y) {
    float a = x[0]-y[0];
    float b = x[1]-y[1];
    return a*a+b*b;
}

int mindist_test(vector x, vector *y, int n) {
    float best = dist(x, y[0]);
    int bestidx = 0;
    for(int i = 1; i &lt; n; i++) {
        float d = dist(x, y[i]);
        if(best &gt; d) {
            best = d;
            bestidx = i;
        }
    }
    return bestidx;
}

typedef struct cell cell;
struct cell {
    vector v;
    cell *next;
};

cell *grid[W][W];
cell cells[N];

float minx, miny, scalex, scaley, scale, scale2;

void build() {
    // Find bounds
    minx = b[0][0];
    miny = b[0][1];
    float maxx = b[0][0];
    float maxy = b[0][1];
    for(int i = 1; i &lt; N; i++) {
        if(minx &gt; b[i][0])
            minx = b[i][0];
        if(maxx &lt; b[i][0])
            maxx = b[i][0];
        if(miny &gt; b[i][1])
            miny = b[i][1];
        if(maxy &lt; b[i][1])
            maxy = b[i][1];
    }
    scalex = (maxx - minx + .00001) / W;
    scaley = (maxy - miny + .00001) / W;
    scale = scalex &lt; scaley ? scalex : scaley;
    scale2 = scale * scale;

    // Place B vectors in grid
    for(int i = 0; i &lt; N; i++) {
        cells[i].v[0] = b[i][0];
        cells[i].v[1] = b[i][1];
        int x = int((b[i][0] - minx) / scalex);
        int y = int((b[i][1] - miny) / scaley);
        cells[i].next = grid[x][y];
        grid[x][y] = &amp;cells[i];
    }
}

int find(vector v) {
    // Start checking at grid[x][y]
    int x = int((v[0] - minx) / scalex);
    int y = int((v[1] - miny) / scaley);

    float mindist = 1e30;
    cell *mincell;

    // Look in grid cells within radius r
    for(int r = 1; r &lt; W; r++) {
        for(int i = x-r; i &lt;= x+r; i++) {
            if(i &lt; 0 || i &gt;= W) continue;
            for(int j = y-r; j &lt;= y+r; j++) {
                if(j &lt; 0 || j &gt;= W) continue;
                for(cell *c = grid[i][j]; c != NULL; c = c-&gt;next) {
                    float d = dist(v, c-&gt;v);
                    if(mindist &gt; d) {
                        mindist = d;
                        mincell = c;
                    }
                }
            }
        }
        // If we've found a point close enough, stop.
        // The next iterations' points are guaranteed to be at least distance r*scale away.
        if(mindist &lt; (r*r)*scale2)
            break;
    }
    return mincell - cells;
}

void thr(int t) {
    int base = N/T*t;
    for(int i = 0; i &lt; N/T; i++) {
        min[base+i] = find(a[base+i]);
    }
}

int main(int argc, char *argv[]) {
    // Initialize
    for(int i = 0; i &lt; N; i++) {
        a[i][0] = rand() / (RAND_MAX/1000000.);
        a[i][1] = rand() / (RAND_MAX/1000000.);
        b[i][0] = rand() / (RAND_MAX/1000000.);
        b[i][1] = rand() / (RAND_MAX/1000000.);
    }

    // Run &amp; time
    double start = mytime();
    build();
    std::thread* threads[T];
    for(int t = 0; t &lt; T; t++) {
        threads[t] = new std::thread(thr, t);
    }
    for(int t = 0; t &lt; T; t++) {
        threads[t]-&gt;join();
    }
    double stop = mytime();
    printf("time: %f\n", stop - start);

    // Check results
    for(int i = 0; i &lt; N; i++) {
        int j = min[i];
        int k = mindist_test(a[i], b, N);
        if(dist(a[i], b[j]) != dist(a[i], b[k])) {
            printf("BAD %d: %d %f %f\n", i, j, dist(a[i], b[j]), dist(a[i], b[k]));
        }
    }
}


Compile with

g++ -O3 -std=c++11 -stdlib=libc++ grid.cc -o grid