Вычислительная матрица парных углов между двумя массивами точек

Question

У меня есть два вектора очков, x и y, в форме (n, p) и (m, p) соответственно. В качестве примера:

x = np.array([[ 0.  , -0.16341, 0.98656], 
       [-0.05937, -0.25205, 0.96589], 
       [ 0.05937, -0.25205, 0.96589], 
       [-0.11608, -0.33488, 0.93508], 
       [ 0.  , -0.33416, 0.94252]]) 
y = np.array([[ 0.  , -0.36836, 0.92968], 
       [-0.12103, -0.54558, 0.82928], 
       [ 0.12103, -0.54558, 0.82928]]) 

Я хочу, чтобы вычислить (n, m) -sized матрицы, содержащих углы между двумя точками, а-ля this вопросом. То есть, Векторизованная версия:

theta = np.array(
      [ np.arccos(np.dot(i, j)/(la.norm(i) * la.norm(j))) 
       for i in x for j in y ] 
     ).reshape((n, m)) 

Примечание: n и m могут быть порядка ~ 10000 каждые.

Answer 1

Есть несколько способов сделать это:

import numpy.linalg as la 
from scipy.spatial import distance as dist 

# Manually 
def method0(x, y): 
    dotprod_mat = np.dot(x, y.T) 
    costheta = dotprod_mat/la.norm(x, axis=1)[:, np.newaxis] 
    costheta /= la.norm(y, axis=1) 
    return np.arccos(costheta) 

# Using einsum 
def method1(x, y): 
    dotprod_mat = np.einsum('ij,kj->ik', x, y) 
    costheta = dotprod_mat/la.norm(x, axis=1)[:, np.newaxis] 
    costheta /= la.norm(y, axis=1) 
    return np.arccos(costheta) 

# Using scipy.spatial.cdist (one-liner) 
def method2(x, y): 
    costheta = 1 - dist.cdist(x, y, 'cosine') 
    return np.arccos(costheta) 

# Realize that your arrays `x` and `y` are already normalized, meaning you can 
# optimize method1 even more 
def method3(x, y): 
    costheta = np.einsum('ij,kj->ik', x, y) # Directly gives costheta, since 
              # ||x|| = ||y|| = 1 
    return np.arccos(costheta) 

результаты синхронизации для (п, т) = (1212, 252):

>>> %timeit theta = method0(x, y) 
100 loops, best of 3: 11.1 ms per loop 
>>> %timeit theta = method1(x, y) 
100 loops, best of 3: 10.8 ms per loop 
>>> %timeit theta = method2(x, y) 
100 loops, best of 3: 12.3 ms per loop 
>>> %timeit theta = method3(x, y) 
100 loops, best of 3: 9.42 ms per loop 

Разница в сроках уменьшает как количество элементы увеличиваются. Для (п, т) = (6252, 1212):

>>> %timeit -n10 theta = method0(x, y) 
10 loops, best of 3: 365 ms per loop 
>>> %timeit -n10 theta = method1(x, y) 
10 loops, best of 3: 358 ms per loop 
>>> %timeit -n10 theta = method2(x, y) 
10 loops, best of 3: 384 ms per loop 
>>> %timeit -n10 theta = method3(x, y) 
10 loops, best of 3: 314 ms per loop 

Однако, если оставить на np.arccos шаг, то есть, предположим, что вы могли бы обойтись только costheta, и не потребность сама theta, затем:

>>> %timeit costheta = np.einsum('ij,kj->ik', x, y) 
10 loops, best of 3: 61.3 ms per loop 
>>> %timeit costheta = 1 - dist.cdist(x, y, 'cosine') 
10 loops, best of 3: 124 ms per loop 
>>> %timeit costheta = dist.cdist(x, y, 'cosine') 
10 loops, best of 3: 112 ms per loop 

Это для случая (6252, 1212). Так что фактически np.arccos занимает 80% времени. В этом случае я считаю, что np.einsum является много быстрее, чем dist.cdist. Поэтому вы определенно хотите использовать einsum.

Резюме: Результаты по theta во многом схожи, но np.einsum является самым быстрым для меня, особенно когда вы не чужеродно вычисления нормы. Старайтесь не вычислять theta и работать с costheta.

Примечание: Важный момент, который я не упомянул о том, что конечность с плавающей точкой точности может вызвать np.arccos дать nan значения. method[0:3] работал для значений x и y, которые не были правильно нормализованы, естественно. Но method3 дал несколько nan s. Я исправил это с предварительной нормализацией, которая, естественно, уничтожает любое усиление при использовании method3, если вам не нужно много раз делать это вычисление для небольшого набора предварительно нормированных матриц (по какой-либо причине).