Как найти попарные различия между строками двух очень больших матриц с помощью numpy?

Учитывая две матрицы, я хочу вычислить попарные различия между всеми строками. Каждая матрица имеет 1000 строк и 100 столбцов, поэтому они довольно большие. Я пробовал использовать цикл for и чистое вещание, но цикл for, похоже, работает быстрее. Я делаю что-то неправильно? Вот код:

from numpy import *
A = random.randn(1000,100)
B = random.randn(1000,100)

start = time.time()
for a in A:
   sum((a - B)**2,1)
print time.time() - start

# pure broadcasting
start = time.time()
((A[:,newaxis,:] - B)**2).sum(-1)
print time.time() - start

Метод вещания длится примерно на 1 секунду дольше, а для больших матриц - еще дольше. Есть идеи, как ускорить это с помощью numpy?


python numpy matrix array-broadcasting
person Bluegreen17    schedule 24.03.2017    source источник

Ответы (4)


arrow_upward
5
arrow_downward

Вот еще один способ выполнения:

(a-b)^2 = a^2 + b^2 - 2ab

с np.einsum для первых двух условий и dot-product для третий -

import numpy as np

np.einsum('ij,ij->i',A,A)[:,None] + np.einsum('ij,ij->i',B,B) - 2*np.dot(A,B.T)

Тест во время выполнения

Подходы -

def loopy_app(A,B):
    m,n = A.shape[0], B.shape[0]
    out = np.empty((m,n))
    for i,a in enumerate(A):
       out[i] = np.sum((a - B)**2,1)
    return out

def broadcasting_app(A,B):
    return ((A[:,np.newaxis,:] - B)**2).sum(-1)

# @Paul Panzer's soln
def outer_sum_dot_app(A,B):
    return np.add.outer((A*A).sum(axis=-1), (B*B).sum(axis=-1)) - 2*np.dot(A,B.T)

# @Daniel Forsman's soln
def einsum_all_app(A,B):
    return np.einsum('ijk,ijk->ij', A[:,None,:] - B[None,:,:], \
                                        A[:,None,:] - B[None,:,:])

# Proposed in this post
def outer_einsum_dot_app(A,B):
    return np.einsum('ij,ij->i',A,A)[:,None] + np.einsum('ij,ij->i',B,B) - \
                                                            2*np.dot(A,B.T)

Сроки -

In [51]: A = np.random.randn(1000,100)
    ...: B = np.random.randn(1000,100)
    ...: 

In [52]: %timeit loopy_app(A,B)
    ...: %timeit broadcasting_app(A,B)
    ...: %timeit outer_sum_dot_app(A,B)
    ...: %timeit einsum_all_app(A,B)
    ...: %timeit outer_einsum_dot_app(A,B)
    ...: 
10 loops, best of 3: 136 ms per loop
1 loops, best of 3: 302 ms per loop
100 loops, best of 3: 8.51 ms per loop
1 loops, best of 3: 341 ms per loop
100 loops, best of 3: 8.38 ms per loop
person Divakar    schedule 24.03.2017
comment
Ха! Обыграл Дивакар на einsum ответ! Конечно, мне кажется, что это неправильный способ его использования, если вам нужно более быстрое решение. . . - person Daniel F; 24.03.2017
comment
@DanielForsman Тебе просто нужно больше практики, в конце концов ты добьешься цели! :) - person Divakar; 24.03.2017
comment
Учитывая, насколько большая часть моего текущего кода полагается на быстрое вычисление декартовых расстояний, этот трюк с алгеброй достаточно полезен, и я не особо возражаю :) - person Daniel F; 24.03.2017

arrow_upward
2
arrow_downward

Вот решение, которое позволяет избежать как цикла, так и больших промежуточных продуктов:

from numpy import *
import time

A = random.randn(1000,100)
B = random.randn(1000,100)

start = time.time()
for a in A:
   sum((a - B)**2,1)
print time.time() - start

# pure broadcasting
start = time.time()
((A[:,newaxis,:] - B)**2).sum(-1)
print time.time() - start

#matmul
start = time.time()
add.outer((A*A).sum(axis=-1), (B*B).sum(axis=-1)) - 2*dot(A,B.T)
print time.time() - start

Печать:

0.546781778336
0.674743175507
0.10723400116
person Paul Panzer    schedule 24.03.2017

arrow_upward
2
arrow_downward

Другая работа для np.einsum

np.einsum('ijk,ijk->ij', A[:,None,:] - B[None,:,:], A[:,None,:] - B[None,:,:])
person Daniel F    schedule 24.03.2017

arrow_upward
0
arrow_downward

Подобно @ paul-panzer, общий способ вычисления попарных разностей массивов произвольной размерности заключается в следующем:

Пусть v будет массивом NumPy размером (n, d):

import numpy as np
v_tiled_across = np.tile(v[:, np.newaxis, :], (1, v.shape[0], 1))
v_tiled_down = np.tile(v[np.newaxis, :, :], (v.shape[0], 1, 1))
result = v_tiled_across - v_tiled_down

Чтобы лучше понять, что происходит, представьте, что каждый d-мерный ряд буквы v приподнят, как флагшток, и скопирован поперек и вниз. Теперь, когда вы выполняете покомпонентное вычитание, вы получаете каждую попарную комбинацию.

__

Также есть scipy.spatial.distance.pdist, который вычисляет метрику попарно.

from scipy.spatial.distance import pdist, squareform
pairwise_L2_norms = squareform(pdist(v))
person maurice    schedule 30.04.2019