Эффективное вычисление евклидова расстояния в python для миллионов строк

Question

Я пытаюсь найти евклидово расстояние между элементами двух наборов данных. У каждого есть миллионы элементов. После расчета евклидова расстояния мне нужно самое близкое совпадение. Учитывая количество элементов, потребуется несколько дней, чтобы закончить

Ниже код, который я пытаюсь. Я также попытался использовать расстояние от scipy.spatial. Но даже это занимает вечно

from sklearn.metrics.pairwise import euclidean_distances
df =pd.DataFrame(euclidean_distances(df1,df2))
df.index =  df1.index
df.columns = df2.index
df['min_distance'] = df.min(axis=1)
df['min_distance_id'] = df.idxmin(axis=1)

Есть ли другой способ получить результат за меньшее время.

Answer 1

Вы смотрели на scipy.spatial.cKDTree?

Вы можете построить эту структуру данных для одного из ваших наборов данных и запросить ее, чтобы получить расстояние для каждой точки во втором наборе данных.

KDTree = scipy.spatial.cKDTree(df1)
distances, indexes = KDTree.query(df2, n_jobs=-1)

Я установил здесь n_jobs=-1, чтобы использовать все доступные процессоры.

Answer 2

Я написал это решение для 2D списков точек, используя numpy. Он быстро найдет ближайшую пару точек между двумя массивами точек. Я попробовал это с двумя списками по 10 миллионов баллов каждый и получил ответ примерно за 4 минуты. С 2 миллионами очков на каждой стороне это заняло всего 42 секунды. Я не знаю, будет ли это достаточно хорошо для ваших нужд, но это определенно быстрее, чем "дни". Это также дает хорошую производительность для больших размеров, если вам это нужно.

def closest(A,B):

    def bruteForce(A,B):
        d = None
        swap = A.shape[0] > B.shape[0]
        if swap: A,B = B,A
        for pA in A:
            daB  = np.sum((pA-B)**2,axis=1)
            iMin = np.argmin(daB)
            if d is None or daB[iMin] < d:
                a,b = pA,B[iMin]
                d   = sum((a-b)**2)
        if swap: a,b = b,a
        return a,b,sqrt(d)

    # small sizes are faster using brute force
    if A.shape[0] * B.shape[0] < 1000000 \
    or A.shape[0] < 20 or B.shape[0] < 20:
        return bruteForce(A,B)

    # find center position
    midA  = np.sum(A,axis=0)/A.shape[0]
    midB  = np.sum(B,axis=0)/B.shape[0]
    midAB = (midA+midB)/2

    # closest A to center position
    A2midAB  = np.sum((A-midAB)**2,axis=1)
    iA       = np.argmin(A2midAB)    
    pA       = A[iA]

    # closest B to pA
    B2pA     = np.sum((B-pA)**2,axis=1)
    iB       = np.argmin(B2pA)
    pB       = B[iB]
    dAB      = sqrt(sum((pA-pB)**2))

    # distance of zero is best solution, return immediately
    if dAB == 0: return pA,pB,dAB

    # slope of ptA-ptB segment
    if pA[0] == pB[0]: p,m = 0,1 
    else:              p,m = 1,(pB[1]-pA[1])/(pB[0]-pA[0])

    # perpendicular line intersections with x axis from each point
    xA = m*A[:,1] + p*A[:,0] 
    xB = m*B[:,1] + p*B[:,0]

    # baselines for ptA and ptB
    baseA = xA[iA]
    baseB = xB[iB]
    rightSide = (baseB > baseA) 

    # partitions
    ArightOfA = (xA > baseA) == rightSide
    BrightOfA = (xB > baseA) == rightSide
    AleftOfB  = (xA > baseB) != rightSide
    BleftOfB  = (xB > baseB) != rightSide

    # include pB and exclude pA (we already know its closest point in B)
    ArightOfA[iA] = False
    AleftOfB[iA]  = False
    BleftOfB[iB]  = True
    BrightOfA[iB] = True

    # recurse left side
    if np.any(AleftOfB) and np.any(BleftOfB):
        lA,lB,lD = closest(A[AleftOfB],B[BleftOfB])
        if lD < dAB: pA,pB,dAB = lA,lB,lD

    # resurse right side
    if np.any(ArightOfA) and np.any(BrightOfA):
        rA,rB,rD = closest(A[ArightOfA],B[BrightOfA])
        if rD < dAB: pA,pB,dAB = rA,rB,rD

    return pA,pB,dAB

Протестировано с использованием двух случайных наборов двумерных точек по 10 миллионов точек каждый:

dimCount = 2
ACount   = 10000000
ASpread  = ACount
BCount   = ACount-1
BSpread  = BCount
A = np.random.random((ACount,dimCount))*ASpread-ASpread/2
B = np.random.random((BCount,dimCount))*BSpread-BSpread/2

a,b,d = closest(A,B)
print("closest points:",a,b,"distance:",d)

# closest points: [-4422004.2963273   2783038.35968559] [-4422004.76974851  2783038.61468366] distance: 0.5377282447465505

То, как это работает, заключается в разделении точек A и B на основе стратегически выбранной пары (pA, pB). Линия между pA и pB служит разделением для точек двух списков. Каждая сторона этого разбиения затем используется рекурсивно для поиска других (более близких) пар точек.

Графически это соответствует разделу на основе перпендикулярных линий сегмента pA-pB:

Стратегия выбора pA и pB состоит в том, чтобы найти приблизительный центр двух групп точек и выбрать точку (pA) из списка A, который находится близко к этому центру. Затем выберите ближайшую точку к pA в списке B. Это гарантирует, что между двумя перпендикулярными линиями, которые ближе к pA или pB в другом списке, нет точек.

Точки A и B, которые находятся на противоположных сторонах перпендикулярных линий, обязательно находятся дальше друг от друга, чем pA-pB, поэтому их можно выделить в двух подсписках и обработать отдельно.

Это позволяет использовать подход «разделяй и властвуй», который значительно сокращает количество расстояний между точками для сравнения.

В моих тестах (со случайно распределенными точками) производительность казалась линейной пропорционально общему количеству точек в A и B. Я попытался исказить распределение, создав небольшие кластеры точек далеко друг от друга (так, чтобы ни одна точка не была на самом деле вблизи приблизительного центра) и производительность была еще линейной. Я не уверен, что существуют какие-либо «наихудшие» распределения точек, которые приведут к падению производительности (я пока не нашел)