Paarweise Vdot mit Numpy

Question

Ich versuche, die paarweise np.vdot eines komplexen 2D-Array x mit sich zu berechnen. Also das Verhalten, das ich will, ist:

X = np.empty((x.shape[0], x.shape[0]), dtype='complex128') 
for i in range(x.shape[0]): 
    for j in range(x.shape[0]): 
     X[i, j] = np.vdot(x[i], x[j]) 

Gibt es eine Möglichkeit, dies ohne die expliziten Schleifen zu tun? Ich versuchte, pairwise_kernel von sklearn zu verwenden, aber es nimmt an, dass die Eingabearrays reelle Zahlen sind. Ich habe auch versucht zu senden, aber vdot flacht seine Eingaben ab.

Answer 1

X = np.einsum('ik,jk->ij', np.conj(x), x) 

entspricht

X = np.empty((x.shape[0], x.shape[0]), dtype='complex128') 
for i in range(x.shape[0]): 
    for j in range(x.shape[0]): 
     X[i, j] = np.vdot(x[i], x[j]) 

np.einsum nimmt eine Summe von Produkten. Der tiefgestellte Index 'ik,jk->ij'np.einsum mitteilt, dass das zweite Argument, np.conj(x) ist ein Array mit Subskripte ik und dem dritten Argumente, x hat Subskripte jk. Somit wird das Produkt np.conj(x)[i,k]*x[j,k] für alle i, j, k berechnet. Die Summe wird über den wiederholten Index k genommen, und da die die i und j übriglässt, werden sie zu den Indizes der resultierenden Anordnung.

---

Zum Beispiel

import numpy as np 

N, M = 10, 20 
a = np.random.random((N,M)) 
b = np.random.random((N,M)) 
x = a + b*1j 

def orig(x): 
    X = np.empty((x.shape[0], x.shape[0]), dtype='complex128') 
    for i in range(x.shape[0]): 
     for j in range(x.shape[0]): 
      X[i, j] = np.vdot(x[i], x[j]) 
    return X 

def alt(x): 
    return np.einsum('ik,jk->ij', np.conj(x), x) 

assert np.allclose(orig(x), alt(x)) 

---

In [307]: %timeit orig(x) 
10000 loops, best of 3: 143 µs per loop 

In [308]: %timeit alt(x) 
100000 loops, best of 3: 8.63 µs per loop 

Answer 2

Um die np.vdot auf alle Zeilen erstrecken, Sie np.tensordot verwenden können, und ich bin zu leihen die konjugierte Idee gerade aus @unutbu's solution, wie so -

np.tensordot(np.conj(x),x,axes=(1,1)) 

Grundsätzlich geben wir mit np.tensordot die zu reduzierenden Achsen an, die in diesem Fall die letzte Achse für die konjugierte Version x und das Array selbst sind, wenn sie auf diese beiden angewendet werden.

Runtime Test -

Lassen Sie uns Zeit @unutbu's solution with np.einsum und die vorgeschlagene Lösung in diesem Beitrag -

In [27]: import numpy as np # From @unutbu's` solution again 
    ...: 
    ...: N, M = 1000, 1000 
    ...: a = np.random.random((N,M)) 
    ...: b = np.random.random((N,M)) 
    ...: x = a + b*1j 
    ...: 

In [28]: %timeit np.einsum('ik,jk->ij', np.conj(x), x) # @unutbu's` solution 
1 loops, best of 3: 4.45 s per loop 

In [29]: %timeit np.tensordot(np.conj(x),x,axes=(1,1)) 
1 loops, best of 3: 3.76 s per loop