effizient parafac/CP Produkt in numpy

Question

Diese Frage konzentriert sich auf numpy.

Ich habe eine Reihe von Matrizen, die alle die gleiche Anzahl von Spalten teilen und unterschiedliche Anzahl von Zeilen haben. Nennen wir sie A, B, C, D usw. und lassen ihre Dimensionen IaxK IbxK, IcxK, usw. sein. ib, ic, id, dh, ...) = \ sum_k A (ia, k) B (ib, k) C (ic, k) ...

Also wenn ich zwei Faktoren habe, lande ich mit einfachem Matrixprodukt.

Natürlich habe ich diese „von Hand“ durch äußere Produkte berechnen kann, so etwas wie:

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 

Dennoch würde ich etwas viel rechentechnisch effizienter mögen, wie unter Berufung numpy Funktionen auf gebautet, aber ich kann nicht finden Relevante Funktionen, kann mir jemand helfen?

Cheers, dank

Answer 1

Vielen Dank sehr viel für Ihre Antworten, habe ich den Tag auf dieser und ich verbrachte schließlich die Lösung gefunden, so dass ich es hier posten für den Datensatz

Diese Lösung erfordert numpy 1.6 und die Verwendung von einsum macht, die Magie leistungsstarke Voodoo ist

im Grunde, wenn Sie hatte Faktor = [A, B, C, D] mit A, B, C und D Matrizen mit die gleiche Anzahl von Spalten, dann würden Sie das PARAFAC Modell berechnen mit:

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

Also, eine Zeile!

Im allgemeinen Fall, habe ich am Ende mit auf den Punkt:

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

im Auge zu haben, dass die äußeree Produkt Kronecker-Produkt in der Verkleidung ist Ihr Problem durch diese einfachen Funktionen gelöst werden soll:

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 nimmt Liste der reellen Zahlen und die Liste der Matrizen.

Nach Kommentar von Jaime bearbeitet.

Answer 3

Ok, so funktioniert das folgende. Zunächst wird ein arbeitet Beispiel dafür, was los ist ...

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) 

Und eine volle Funktion

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)