Wie finden Sie die paarweisen Unterschiede zwischen Zeilen von zwei sehr großen Matrizen mit numpy?

Question

Gegeben zwei Matrizen, ich möchte die paarweisen Unterschiede zwischen allen Zeilen berechnen. Jede Matrix hat 1000 Zeilen und 100 Spalten, so dass sie ziemlich groß sind. Ich habe versucht, eine for-Schleife und reine Rundfunkübertragung zu verwenden, aber die for-Schleife scheint schneller zu arbeiten. Mache ich etwas falsch? Hier ist der Code:

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 

Die Broadcasting-Methode dauert etwa 1 Sekunde länger und es ist noch länger für große Matrizen. Irgendeine Idee, wie man das mit numpy rein beschleunigt?

Answer 1

Hier ist eine andere Art und Weise durchzuführen :

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

mit np.einsum für die ersten zwei Terme und dot-product für den dritten -

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) 

Runtime Test

Approaches -

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) 

Zeiten -

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

Hier ist eine Lösung, die sowohl die Schleife und die großen Zwischenprodukte vermeidet:

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 

Drucke:

0.546781778336 
0.674743175507 
0.10723400116 

Answer 3

Ein weiterer Job für np.einsum

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