Czy numpy zapewnia uogólniony wewnętrzny produkt?

Question

Większość języków zorientowanych na tablice, takich jak APL lub J, ma pewną formę uogólnionego produktu wewnętrznego, który może działać jak standardowe mnożenie macierzy, ale obsługuje arbitralne operacje w miejsce standardowych. Na przykład w J +/ . * jest standardem mnożenia, a następnie sumą, ale możesz też np. <./ . +, aby uzyskać operację dodawania, a następnie min (powiedz na temat przyrostowej aktualizacji długości najkrótszej ścieżki przez wykres).

W powolne i 2D-tylko Pythonie, byłoby to coś jak:

import numpy as np 

def general_inner(f, g, x, y): 
    return np.array([[f(g(x1, y1)) for y1 in y.T] for x1 in x]) 

x = np.arange(1, 5, dtype="float").reshape((2, 2)) 
y = np.array([[0.9], [0.1]]) 
assert(x.dot(y) == general_inner(np.sum, np.multiply, x, y)) 

Czy numpy zapewnić wszystko, aby bezpośrednio obsługiwać ten wzór?

Answer 1

Możesz dostać się tam z krojeniem. Możemy przekształcić te dwa argumenty tak, aby operacja otrzymała rozgłaszanie zamiast wykonywanej elementwise, a następnie wykonać operację redukcji wzdłuż niechcianej osi.

import numpy 

a = numpy.array([[1, 2, 3], 
       [4, 5, 6]]) 
b = numpy.array([[7, 8], 
       [9, 10], 
       [11, 12]]) 

# Ordinary matrix multiplication 
print(a @ b) 
# Matrix multiplication using broadcasting 
print(numpy.sum(a[:,:,numpy.newaxis] * b[numpy.newaxis,:,:], axis=1)) 
# Our "generalized" version 
print(numpy.min(a[:,:,numpy.newaxis] + b[numpy.newaxis,:,:], axis=1)) 

wahałbym się, że jest to „produkt uogólniony wewnętrzna” call, ponieważ produkty wewnętrzne mają specyficzną strukturę matematyczną, że ta nowa wersja brakuje.

Sposób ten w zasadzie działa to, że każdy numpy.newaxis ma długość 1 i dostaje nadawanie, więc:

a[:,:,numpy.newaxis] * b[numpy.newaxis,:,:] 

daje nam:

result[i,j,k] = a[i,j] * b[j,k] 

Albo jeśli to pomaga zrozumieć (znajdę nadawanie czasami nieco mylące),

aa = a[:,:,numpy.newaxis] 
bb = b[numpy.newaxis,:,:] 
result[i,j,k] = aa[i,j,0] * bb[0,j,k] 

Answer 2

Nie tak wolno numpy odpowiednik to g(x,y.T), z wykorzystaniem nadawania, a następnie f(..., axis=1).

In [136]: general_inner(np.sum, np.multiply, x, y) 
Out[136]: 
array([[ 1.1], 
     [ 3.1]]) 
In [137]: np.multiply(x,y.T) 
Out[137]: 
array([[ 0.9, 0.2], 
     [ 2.7, 0.4]]) 
In [138]: np.sum(np.multiply(x,y.T),axis=1) 
Out[138]: array([ 1.1, 3.1]) 

podobnie do maksymalnej kwoty:

In [145]: general_inner(np.max, np.add, x, y) 
Out[145]: 
array([[ 2.1], 
     [ 4.1]]) 
In [146]: np.max(np.add(x,y.T), axis=1) 
Out[146]: array([ 2.1, 4.1]) 

Istnieje potencjalne zamieszanie w które np.add, np.multiply, np.maximum są ufunc, natomiast np.sum, np.prod, np.max nie są, ale warto parametru axis, i keepdims. (Edit: np.add.reduce jest ufunc odpowiednik np.sum.)

In [152]: np.max(np.add(x,y.T), axis=1, keepdims=True) 
Out[152]: 
array([[ 2.1], 
     [ 4.1]]) 

Nie był stary prośba akcesorium do wdrożenia tego rodzaju rzeczy w ciągu np.einsum. To implementuje obliczenia sum of products z wysokim stopniem kontroli nad indeksowaniem. Więc koncepcyjnie może wykonać max of sums z taką samą kontrolą indeksowania. Ale o ile wiem, nikt nie próbował go wdrożyć.

Ten uogólniony produkt wewnętrzny był uroczą cechą APL (to był mój podstawowy język kilka dekad temu). Ale najwyraźniej nie jest tak przydatny, że migrował z tej rodziny języków. MATLAB nie ma czegoś takiego.

Czy jest coś, co APL & J może zrobić z tym konstruktem, czego nie można dokonać w numpy z rodzajem nadawania, które demonstrowaliśmy?

---

Z bardziej ogólnych x i y kształtach, muszę dodać newaxis jak podano w drugiej odpowiedzi

In [176]: x = np.arange(3*4).reshape(4,3) 
In [177]: y = np.arange(3*2).reshape(3,2) 
In [178]: np.sum(np.multiply(x[...,None],y[None,...]),axis=1) 
Out[178]: 
array([[10, 13], 
     [28, 40], 
     [46, 67], 
     [64, 94]]) 
In [179]: np.max(np.add(x[...,None],y[None,...]),axis=1) 
Out[179]: 
array([[ 6, 7], 
     [ 9, 10], 
     [12, 13], 
     [15, 16]]) 

uogólniając do 3D, z wykorzystaniem matmul ideę ostatni dim/2nd ostatni z matmul:

In [195]: x = np.arange(2*4*5).reshape(2,4,5) 
In [196]: y = np.arange(2*5*3).reshape(2,5,3) 
In [197]: np.einsum('ijk,ikm->ijm', x, y).shape 
Out[197]: (2, 4, 3) 
In [203]: np.add.reduce(np.multiply(x[...,None], y[...,None,:,:]), axis=-2).shape 
Out[203]: (2, 4, 3) 

# shapes broadcast: (2,4,5,n) * (2,n,5,3) => (2,4,5,3); sum on the 5 

Więc gdy numpy (i MATLAB) nie ma specjalnego Składnia taka jak APL, pomysł operacji rozszerzania (na wyjściu), po której następuje redukcja, jest dość powszechny.

testowania innych ufunc:

In [205]: np.maximum.reduce(np.add(x[...,None], y[...,None,:,:]), axis=-2).shape 
Out[205]: (2, 4, 3) 
In [208]: np.logical_or.reduce(np.greater(x[...,None], y[...,None,:,:]), axis=-2).shape 
Out[208]: (2, 4, 3)