Расчет расстояний между векторами двух numpy массивов

Question

У меня есть два numpy массива R с размерами S x F и W с размерами N x M x F . Получение бетона позволяет присвоить следующие значения N = 5, M = 7, F = 3, S = 4

Массив R содержит коллекции образцов S = 4 с F = 3 функциями , Каждая строка представляет образцы, а каждая строка - объект. Поэтому R[0] - это первый образец, R[1] - второй и продолжается. Каждая запись R[i-th] содержит F элементов, например, 1023 *.

Вот небольшой фрагмент для инициализации всех этих значений с учетом MWE

import numpy as np

# Size of Map (rows, columns)
N, M = 5, 7

# Number of features
F = 3

# Sample size
S = 4

np.random.seed(13)
R = np.random.randint(0, 10, size=(S, F))
W = np.random.randint(-4, 5, size=(N, M, F))

Мы также можем видеть заданную «линию глубины» numpy массив W , как вектор, также имеющий тот же размер, что и каждая строка массива R (это легко заметить, глядя на размер последнего измерения обоих массивов). С этим я могу получить доступ к W[2, 3] и получить np.array([ 2, 2, -1 ]) (значения здесь только примеры).

Я создал простую функцию для вычисления расстояния данного вектора r до каждого «линия глубины» матрицы W и возврат положения ближайшего элемента W линии глубины в r

def nearest_vector_matrix_naive(r, W):
    delta = np.zeros((N,M), dtype=int)
    for i in range(N):
        for j in range(M):
            norm = 0
            for k in range(F):
                norm += (r[k] - W[i,j,k])**2
            delta[i,j] = norm
            norm = 0
    win_idx = np.unravel_index(np.argmin(delta, axis=None), delta.shape)
    return win_idx

Конечно, это очень наивный подход, который я мог бы дополнительно оптимизировать для приведенного ниже кода, получив ОГРОМНОЕ повышение производительности.

def nearest_vector_matrix(r, W):
    delta = np.sum((W[:,:] - r)**2, axis=2)
    return np.unravel_index(np.argmin(delta, axis=None), delta.shape)

Я могу использовать эту функцию просто как

nearest_idx = nearest_vector_matrix(R[0], W)
# Returns the nearest vector in W to R[0]
W[nearest_idx]

Поскольку у меня есть массив R с набором выборок, я использую следующий фрагмент для вычисления ближайших векторов к массиву выборок:

def nearest_samples_matrix(R, W):
    DELTA = np.zeros((R.shape[0],2))
    for idx, r in enumerate(R):
        delta = np.sum((W[:,:] - r)**2, axis=2)
        DELTA[idx] = np.unravel_index(np.argmin(delta, axis=None), delta.shape)
    return DELTA

Эта функция возвращает массив с S строками ( S - количество выборок) двумерных индексов. То есть DELTA имеет (S, 2) форму (всегда).

Я хотел бы знать, как я могу заменить for l oop (например, для трансляции) внутри nearest_samples_matrix для улучшения кода производительность исполнения еще дальше?

Я не мог понять, как это сделать. (кроме того, я смог сделать это в первом случае)

Answer 1

Лучшее решение зависит от входного размера массивов

Для задач меньшего размера dim <20 или менее, <a href="https://docs.scipy.org/doc/scipy/reference/generated/scipy.spatial.cKDTree.html" rel="nofollow noreferrer"> подход kdtree обычно является способом к go. Есть довольно много ответов относительно этой темы c, например. one Я написал несколько недель: go.

Если размер проблемы слишком велик, вы можете перейти к алгоритмам перебора. Оба следующих алгоритма намного быстрее, чем ваш оптимизированный подход, но при больших входных размерах и задачах с малыми размерами гораздо медленнее, чем подход kdtree O (log (n)) вместо O (n ^ 2).

Грубая сила 1

В следующем примере используется алгоритм, описанный здесь . Это очень быстро для задач большого размера, потому что большая часть вычислений выполняется в высоко оптимизированном алгоритме матрично-матричного умножения. Недостатком является высокий уровень использования памяти (все расстояния рассчитываются за один вызов функции) и проблемы с точностью из-за более подверженного ошибкам метода вычисления.

import numpy as np
from sklearn.metrics.pairwise import euclidean_distances

def nearest_samples_matrix_2(R,W):
    R_Temp=R
    W_Temp=W.reshape(-1,W.shape[2])
    dist=euclidean_distances(R_Temp, W_Temp)
    ind_1,ind_2=np.unravel_index(np.argmin(dist,axis=1),shape=(W.shape[0],W.shape[1]))
    return np.vstack((ind_1,ind_2)).T

Грубая сила 2

Это очень похоже на ваш наивный подход, но использует JIT-компилятор (Numba) для получения хорошей производительности. Временные массивы не нужны, и точность должна быть хорошей (до тех пор, пока не произойдет переполнение). Есть место для дальнейшей оптимизации (l oop tiling) при больших входных размерах.

import numpy as np
import numba as nb

#parallelization is only beneficial on larger input data
@nb.njit(fastmath=True,parallel=True,cache=True)
def nearest_samples_matrix_3(r, W):
    ind_i=0
    ind_j=0
    out=np.empty((r.shape[0],2),dtype=np.int64)
    for x in nb.prange(r.shape[0]):
        delta=0
        for k in range(W.shape[2]):
            delta += (r[x,k] - W[0,0,k])**2

        for i in range(W.shape[0]):
            for j in range(W.shape[1]):
                norm = 0
                for k in range(W.shape[2]):
                    norm += (r[x,k] - W[i,j,k])**2
                if norm < delta:
                    delta=norm
                    ind_i=i
                    ind_j=j
        out[x,0]=ind_i
        out[x,1]=ind_j
    return out

Сроки

#small Arrays
N, M = 100, 200
F = 30
S = 50
R = np.random.randint(0, 10, size=(S, F))
W = np.random.randint(-4, 5, size=(N, M, F))

#your function
%timeit nearest_samples_matrix(R,W)
#268 ms ± 2.94 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
%timeit nearest_samples_matrix_2(R,W)
#5.62 ms ± 22.6 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
%timeit nearest_samples_matrix_3(R,W)
#3.68 ms ± 1.01 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

#larger arrays
N, M = 1_000, 2_000
F = 50
S = 100
R = np.random.randint(0, 10, size=(S, F))
W = np.random.randint(-4, 5, size=(N, M, F))

#%timeit nearest_samples_matrix_1(R,W)
#too slow
%timeit nearest_samples_matrix_2(R,W)
#2.76 s ± 17.2 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
%timeit nearest_samples_matrix_3(R,W)
#1.42 s ± 402 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)