nietrywialnym sumy produktów zewnętrznych bez tymczasowych w numpy

Question

rzeczywisty problem życzę rozwiązania jest, biorąc pod uwagę zestaw N wektor jednostkowy i inny zestaw M wektorów obliczyć dla każdej jednostki wektory średnią wartość bezwzględna iloczynu punktowego z każdym z wektorów M. Zasadniczo polega to na obliczeniu zewnętrznego produktu obu matryc oraz zsumowaniu i uśrednieniu z wartością bezwzględną utkniętą pomiędzy.

Dla N i M nie jest zbyt duży to nie jest trudne i istnieje wiele sposobów, aby kontynuować (patrz niżej). Problem polega na tym, że utworzone pliki tymczasowe są ogromne i stanowią praktyczne ograniczenie dla dostarczonego podejścia. Czy te obliczenia można wykonać bez tworzenia tymczasowych? Główna trudność, jaką mam, wynika z obecności wartości bezwzględnej. Czy istnieją ogólne techniki "gwintowania" takich obliczeń?

Jako przykład należy rozważyć następujący kod

N = 7 
M = 5 

# Create the unit vectors, just so we have some examples, 
# this is not meant to be elegant 
phi = np.random.rand(N)*2*np.pi 
ctheta = np.random.rand(N)*2 - 1 
stheta = np.sqrt(1-ctheta**2) 
nhat = np.array([stheta*np.cos(phi), stheta*np.sin(phi), ctheta]).T 

# Create the other vectors 
m = np.random.rand(M,3) 

# Calculate the quantity we desire, here using broadcasting. 
S = np.average(np.abs(np.sum(nhat*m[:,np.newaxis,:], axis=-1)), axis=0) 

to jest wielki S jest tablicą długości N i zawiera żądane wyniki. Niestety w tym procesie stworzyliśmy potencjalnie ogromne tablice. Wynikiem

np.sum(nhat*m[:,np.newaxis,:], axis=-1) 

jest M X N tablicy. Końcowy wynik, oczywiście, ma tylko rozmiar N. Zacznij zwiększać rozmiary N i M, a my szybko napotkamy błąd pamięci.

Jak stwierdzono powyżej, jeśli wartość bezwzględna nie wymaga to można postępować jak następuje teraz za pomocą einsum()

T = np.einsum('ik,jk,j', nhat, m, np.ones(M))/M 

Działa to działa szybko i nawet w przypadku bardzo dużej N i M. Aby rozwiązać konkretny problem, potrzebuję dodać abs(), ale bardziej ogólne rozwiązanie (być może bardziej ogólne ufunc) również byłoby interesujące.

Answer 1

Na podstawie niektórych komentarzy wydaje się, że najlepszym rozwiązaniem jest użycie cythonu. Głupio nigdy nie patrzyłem na użycie cythonu. Okazuje się, że stosunkowo łatwo jest wyprodukować działający kod.

Po kilku poszukiwaniach złożę następujący kod cytonu. Jest to najbardziej ogólny kod, prawdopodobnie nie jest najlepszym sposobem na napisanie kodu, ale prawdopodobnie może być bardziej efektywny. Mimo to jest to tylko około 25% wolniej niż kod einsum() w oryginalnym pytaniu, więc nie jest tak źle! Został napisany tak, by działał jawnie z tablicami utworzonymi tak, jak zrobiono w pierwotnym pytaniu (stąd założone tryby tablic wejściowych).
Pomimo zastrzeżeń, które stanowią rozsądnie skuteczne rozwiązanie pierwotnego problemu i może służyć jako punkt wyjścia w podobnych sytuacjach.

import numpy as np 
cimport numpy as np 
import cython 
DTYPE = np.float64 
ctypedef np.float64_t DTYPE_t 
cdef inline double d_abs (double a) : return a if a >= 0 else -a 

@cython.boundscheck(False) 
@cython.wraparound(False) 
def process_vectors (np.ndarray[DTYPE_t, ndim=2, mode="fortran"] nhat not None, 
        np.ndarray[DTYPE_t, ndim=2, mode="c"] m not None) : 
    if nhat.shape[1] != m.shape[1] : 
     raise ValueError ("Arrays must contain vectors of the same dimension") 
    cdef Py_ssize_t imax = nhat.shape[0] 
    cdef Py_ssize_t jmax = m.shape[0] 
    cdef Py_ssize_t kmax = nhat.shape[1] # same as m.shape[1] 
    cdef np.ndarray[DTYPE_t, ndim=1] S = np.zeros(imax, dtype=DTYPE) 
    cdef Py_ssize_t i, j, k 
    cdef DTYPE_t val, tmp 
    for i in range(imax) : 
     val = 0 
     for j in range(jmax) : 
      tmp = 0 
      for k in range(kmax) : 
       tmp += nhat[i,k] * m[j,k] 
      val += d_abs(tmp) 
     S[i] = val/jmax 
    return S 

Answer 2

Jest to nieco wolniej, ale nie tworzy dużej macierzy pośredniej.

vals = np.zeros(N) 
for i in xrange(N): 
    u = nhat[i] 
    for v in m: 
     vals[i]+=abs(np.dot(u,v)) 
    vals[i]=vals[i]/M 

edytuj: przesunięty podział przez M poza pętlą for.

edit2: nowy pomysł, zachowując stary dla potomności i odpowiedni komentarz.

m2 = np.average(m,0) 
vals = np.zeros(N) 
for i in xrange(N): 
    u=nhat[i] 
    vals[i]=abs(np.dot(u,m2)) 

To jest szybkie, ale czasami daje różne wartości, pracuję nad tym, ale może może pomóc w średnim czasie.

edytuj 3: Ah, to rzecz absolutnej wartości.hmm

>>> S 
array([ 0.28620962, 0.65337876, 0.37470707, 0.46500913, 0.49579837, 
     0.29348924, 0.27444208, 0.74586928, 0.35789315, 0.3079964 , 
     0.298353 , 0.42571445, 0.32535728, 0.87505053, 0.25547394, 
     0.23964505, 0.44773271, 0.25235646, 0.4722281 , 0.33003338]) 
>>> vals 
array([ 0.2099343 , 0.6532155 , 0.33039334, 0.45366889, 0.48921527, 
     0.20467291, 0.16585856, 0.74586928, 0.31234917, 0.22198642, 
     0.21013519, 0.41422894, 0.26020981, 0.87505053, 0.1199069 , 
     0.06542492, 0.44145805, 0.08455833, 0.46824704, 0.28483342]) 

time to compute S: 0.000342130661011 seconds 
time to compute vals: 7.29560852051e-05 seconds 

edit 4: Dobrze, jeśli mają wartość głównie pozytywne dla wektorów jednostkowych powinno to działać szybciej, zakładając wektory wm zawsze są pozytywne, jak są one w swoich fikcyjnych danych.

m2 = np.average(m,0) 
vals = np.zeros(N) 
for i in xrange(N): 
    u=nhat[i] 
    if u[0] >= 0 and u[1] >= 0 and u[2] >= 0: 
     vals[i] = abs(np.dot(u,m2)) 
    else: 
     for j in xrange(M): 
      vals[i]+=abs(np.dot(u,m[j])) 
     vals[i]/=M 

Answer 3

Nie sądzę, że istnieje łatwy sposób (poza Cythonem i tym podobne), aby przyspieszyć dokładną operację. Ale możesz zastanowić się, czy naprawdę musisz obliczyć to, co obliczasz. Bo jeśli zamiast średniej wartości bezwzględnych można użyć root mean square, by nadal być jakoś średnio jasności produktów wewnętrznych, ale można dostać go w jednym strzale jak:

rms = np.sqrt(np.einsum('ij,il,kj,kl,k->i', nhat, nhat, m, m, np.ones(M)/M)) 

To jest taka sama jak robi:

rms_2 = np.sqrt(np.average(np.einsum('ij,kj->ik', nhat, m)**2, axis=-1)) 

Tak, to nie jest dokładnie to, co prosiłeś, ale obawiam się, że jest tak blisko, jak dostaniesz z vectorized podejścia. Jeśli zdecydujesz się pójść tą drogą, zobacz, jak dobrze np.einsum wykonuje dla dużych N i M: ma tendencję do trzebieży, gdy przekazał zbyt wiele parametrów i indeksów.