如何找到使用numpy的兩個非常大的矩陣行之間的成對差異？

Question

給定兩個矩陣，我想計算所有行之間的成對差異。每個矩陣有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？

Answer 1

另一個工作下面就來進行另一種方式：

（AB）^ 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 

Answer 2

這裏是避免了兩者的環和大中間體的解決方案：

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 

Answer 3

爲np.einsum

np.einsum('ijk,ijk->ij', A[:,None,:] - B[None,:,:], A[:,None,:] - B[None,:,:])