Fastest array-by-array nearest neighbour search
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