calcolo efficiente del prodotto parafac/CP in numpy

Question

Questa domanda si concentra su numpy.

Ho un set di matrici che condividono tutti lo stesso numero di colonne e hanno un numero diverso di righe. Chiamiamoli A, B, C, D, ecc e lascia che le loro dimensioni siano IaxK IbxK, IcxK, ecc

Quello che voglio è calcolare in modo efficiente il tensore IaxIbxIc ... definito come segue: P (ia, ib, ic, id, ie, ...) = \ sum_k A (ia, k) B (ib, k) C (ic, k) ...

Quindi se ho due fattori, finisco con un prodotto a matrice semplice.

Certo che posso calcolare questo "a mano" attraverso prodotti esterni, qualcosa di simile:

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 

Eppure, vorrei qualcosa di molto più computazionalmente efficiente, come basandosi su builtin funzioni NumPy, ma non riesco a trovare funzioni rilevanti, qualcuno può aiutarmi?

Saluti, grazie

Answer 1

Grazie a tutti voi per le vostre risposte, ho trascorso la giornata su questo e alla fine ho trovato la soluzione, così ho posto qui per la cronaca

Questa soluzione richiede NumPy 1.6 e si avvale di einsum, che è potente voodoo magia

fondamentalmente, se si ha fattore = [A, B, C, D] con A, B, C e D matrici con lo stesso numero di colonne, allora sarebbe calcolare il modello parafac utilizzando:

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

così, una riga!

Nel caso generale, io alla fine con questo:

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

avendo in mente quel prodotto esterno è prodotto Kronecker mascherato il problema dovrebbe essere risolto da questo semplice funzioni:

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 prende elenco di numeri reali e la lista delle matrici.

Modificato dopo il commento di Jaime.

Answer 3

Ok, quindi il seguente funziona. Prima un esempio pratico di ciò che sta succedendo ...

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) 

E una funzione completa

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)