wydajnie obliczanie produktu PARAPAC/CP w numpy

Question

To pytanie skupia się na numpy.

Mam zestaw macierzy, które wszystkie mają taką samą liczbę kolumn i mają różną liczbę wierszy. Nazwijmy je A, B, C, D, itd. Niech ich wymiary będą takie jak IaxK IbxK, IcxK, itp.

To, czego chcę, to wydajnie obliczyć IaxIbxIc ... tensor P zdefiniowany następująco: P (ia, ib, ic, id, ie, ...) = \ sum_k A (ia, k) B (ib, k) C (ic, k) ...

Więc jeśli mam dwa czynniki, kończę z prostym produktem matrycowym.

Oczywiście mogę obliczyć to „ręcznie” przez produktów zewnętrznych, coś jak:

def parafac(factors,components=None): 
     ndims = len(factors) 
     ncomponents = factors[0].shape[1] 
     total_result=array([]) 
     if components is None: 
      components=range(ncomponents) 

     for k in components: 
      #for each component (to save memory) 
      result = array([]) 
      for dim in range(ndims-1,-1,-1): 
       #Augments model with next dimension 
       current_dim_slice=[slice(None,None,None)] 
       current_dim_slice.extend([None]*(ndims-dim-1)) 
       current_dim_slice.append(k) 
       if result.size: 
        result = factors[dim].__getitem__(tuple(current_dim_slice))*result[None,...] 
       else: 
        result = factors[dim].__getitem__(tuple(current_dim_slice)) 
      if total_result.size: 
       total_result+=result 
      else: 
       total_result=result 
     return total_result 

Still, chciałbym coś znacznie bardziej obliczeniowo wydajne, jak opierając się na wbudowane funkcje NumPy, ale nie mogę znaleźć odpowiednie funkcje, czy ktoś może mi pomóc?

Cheers, dzięki

Answer 1

Dziękuję wszystkim bardzo za odpowiedzi, ja spędziłem dzień na to i ostatecznie znaleźć rozwiązanie, więc po to tutaj dla rekordu

to rozwiązanie wymaga numpy 1.6 i wykorzystuje einsum, który jest potężny magia voodoo

Zasadniczo, jeśli miał współczynnik = [A, B, C, D] z A, B, C i D matryc taka sama liczba kolumn, a następnie by obliczyć wzór parafac używany:

import numpy 
P=numpy.einsum('az,bz,cz,dz->abcd',A,B,C,D) 

więc jedna linia!

W ogólnym przypadku, skończę z tym:

def parafac(factors): 
    ndims = len(factors) 
    request='' 
    for temp_dim in range(ndims): 
     request+=string.lowercase[temp_dim]+'z,' 
    request=request[:-1]+'->'+string.lowercase[:ndims] 
    return einsum(request,*factors) 

Answer 2

mając na uwadze, że produkt zewnętrzna jest produktem Kronecker w przebraniu Twój problem powinien być rozwiązany poprzez ten prosty funkcji:

def outer(vectors): 
    shape=[v.shape[0] for v in vectors] 
    return reduce(np.kron, vectors).reshape(shape) 
def cp2Tensor(l,A): 
    terms=[]  
    for r in xrange(A[0].shape[1]): 
     term=l[r]*outer([A[n][:,r] for n in xrange(len(A))]) 
     terms.append(term) 
    return sum(terms) 

cp2Tensor bierze listę liczb rzeczywistych i listy matryc.

Edytowane po komentarzu przez Jaime.

Answer 3

Ok, więc następujące prace. Najpierw pracował na przykład tego, co się dzieje ...

a = np.random.rand(5, 8) 
b = np.random.rand(4, 8) 
c = np.random.rand(3, 8) 
ret = np.ones(5,4,3,8) 
ret *= a.reshape(5,1,1,8) 
ret *= b.reshape(1,4,1,8) 
ret *= c.reshape(1,1,3,8) 
ret = ret.sum(axis=-1) 

i pełną funkcją

def tensor(elems) : 
    cols = elems[0].shape[-1] 
    n_elems = len(elems) 
    ret = np.ones(tuple([j.shape[0] for j in elems] + [cols])) 
    for j,el in enumerate(elems) : 
     ret *= el.reshape((1,) * j + (el.shape[0],) + 
          (1,) * (len(elems) - j - 1) + (cols,)) 
    return ret.sum(axis=-1)