How can I efficiently find the minimum distance values for points in two very large 3D coordinate arrays?

Let's pretend we have two arrays of 3D coordinate points.

A = (1000000, 3) of type float,
B = (100000, 3) of type float

For each coordinate in A, I would like to find the minimum euclidean distance to any coordinate in B. This means that it should calculate the euclidean distance between A[0], for instance, with all of B, and then take the smallest value.

I wrote up code to do this using loops. It works, but due to the size of my arrays it takes more than an hour to complete. The pseudocode looks something like this:

minDistances = np.zeros(A.shape)
for i in range(len(A)):
  queriedPoint = A[i]
  distances = B - queriedPoint
  euclideanDistances = np.linalg.norm(distances, axis=1)
  minDistance = np.min(euclideanDistances)
  minDistances[i] = minDistance

Ideally, I'd like to vectorize this, but in doing so it seems to crash my program from RAM usage. Any ideas or tips for how to do this more efficiently? I was wondering if maybe I could reduce the problem to something easier somehow, or else rethink how to approach it alltogether. Thanks!

The fastest way to do this is probably using scipy.spatial.KDTree. The proposed duplicate recommends using scipy.spatial.distance.cdist, but your arrays are too large that it would consume too much memory.

import numpy as np
from scipy.spatial import KDTree

rng = np.random.default_rng(42)
A = rng.uniform(low=-100, high=100, size=(1_000_000, 3))
B = rng.uniform(low=-100, high=100, size=(100_000, 3))

tree = KDTree(B)
distances = tree.query(A)[0]

I don't know the actual range of values in A and B, so I just used (-100, 100). This takes about 2.2 seconds to run.

Brute force will require NxN distance calculations.

A simpler method is to use "bucketing", i.e. use some special boxes.

Building boxes requires about 4N calculations. For example, first pass determine the max & min X,Y,Z, coords for each array. Then you split the "space" into say 4x4x4=64 boxes for array_1 and other 64 for array_2.

By simple vertices comparisons you can get the two boxes (one for each array) that are closer. Yes, this is brute force, but with small data.
<br>Note: If boxes intersect or there are more than one close pair, then you need a list of candidates, more boxes, but still not the initial big data.

Then run new pass on the arrays. Get only those points that lay in the list of boxes.

Finally you can run a brute-force with that chosen points.

For a pesimal case, where most of boxes of array_1 intersect (or are at same distance) with some boxes of array_2, then just split each box into eight smaller boxes and check again. The worst case may be (but rare) worst that brute force with the arrays.

Yes, so the best way is to divide the problem into subproblems, divide and conquer algorithm so given the pseudocode above, we can try to attempt it by utilizing dict & list comprehension

import numpy as np
A = np.random.rand(1000000,3)
B = np.random.rand(1000000,3)

minDistance = {i: np.linalg.norm(B-A[i],axis=1).min() for i in range(len(A))}
minDistance = [minDistance[i] for i in range(len(A))]

print(minDistance)

by using dict & list comprehension, they help to optimize both time and complexity. hope that helps

Alternatively, you can try the loop-based approach

import numpy as np
from scipy.spatial.distance import cdist

A = np.random.rand(1000000,3)
B = np.random.rand(1000000,3)

minDist = np.zeros(len(A))

for i, coord in enumerate(A):
    distances = cdist(np.expand_dims(coord,axis=0),B)
    minDist[i] = np.min(distances)
print(minDist)

If the two arrays are filled with completely random numbers, then perhaps there’s little to be done. If each array corresponds to for example a vehicle track, then you should look at the Hausdorff distance and the Fréchet metric.