Loop sopra (o vectorize) matrici di lunghezza variabile in Teanò

Question

Ho un elenco di matrici L, in cui ogni elemento è una matrice Mx*n (x è una variabile, n è una costante).

voglio calcolare la somma dei M'*M per tutti gli elementi a L (M' è la trasposizione di M) come il seguente codice Python fa:

for M in L: 
    res += np.dot(M.T, M) 

In realtà voglio implementare questo in Theano (che doesn 't supportano array multidimensionali a lunghezza variabile), e non voglio riempire tutte le matrici con le stesse dimensioni perché ciò sprecerebbe troppo spazio (alcune matrici possono essere molto grandi).

C'è un modo migliore per farlo?

Edit:

L è noto prima della compilazione Theano.

Edit:

ha ricevuto due risposte eccellenti da @DanielRenshaw e @Divakar, emotivamente difficile sceglierne uno da accettare.

Answer 1

È possibile eseguire il rilievo degli array di input lungo il primo asse che è il totale di x. Quindi, ci ritroveremmo con un alto array (X,n), dove X =x1+x2+x3+..... Questo può essere trasposto e il suo prodotto punto con se stesso sarebbe l'output desiderato di forma (n,n). Tutto ciò si ottiene con una soluzione vettoriale pura che sfrutta il potente prodotto dot.Pertanto, l'applicazione sarebbe -

# Concatenate along axis=0 
Lcat = np.concatenate(L,axis=0) 

# Perform dot product of the transposed version with self 
out = Lcat.T.dot(Lcat) 

test runtime e verificare uscita -

In [116]: def vectoized_approach(L): 
    ...: Lcat = np.concatenate(L,axis=0) 
    ...: return Lcat.T.dot(Lcat) 
    ...: 
    ...: def original_app(L): 
    ...: n = L[0].shape[1] 
    ...: res = np.zeros((n,n)) 
    ...: for M in L: 
    ...:  res += np.dot(M.T, M) 
    ...: return res 
    ...: 

In [117]: # Input 
    ...: L = [np.random.rand(np.random.randint(1,9),5)for iter in range(1000)] 

In [118]: np.allclose(vectoized_approach(L),original_app(L)) 
Out[118]: True 

In [119]: %timeit original_app(L) 
100 loops, best of 3: 3.84 ms per loop 

In [120]: %timeit vectoized_approach(L) 
1000 loops, best of 3: 632 µs per loop 

Answer 2

Dato che il numero di matrici è noto prima che la compilazione Theano debba avvenire, è sufficiente utilizzare elenchi Python regolari delle matrici di Theano.

Ecco un esempio completo che mostra la differenza tra le versioni numpy e Theano.

Questo codice è stato aggiornato per includere un confronto con l'approccio vettorializzato di @Divakar che offre prestazioni migliori. Sono possibili due approcci vettorizzati per Theano, uno in cui Theano esegue la concatenazione, e uno in cui numpy fa la concatenazione il cui risultato viene poi passato a Theano.

import timeit 
import numpy as np 
import theano 
import theano.tensor as tt 


def compile_theano_version1(number_of_matrices, n, dtype): 
    assert number_of_matrices > 0 
    assert n > 0 
    L = [tt.matrix() for _ in xrange(number_of_matrices)] 
    res = tt.zeros(n, dtype=dtype) 
    for M in L: 
     res += tt.dot(M.T, M) 
    return theano.function(L, res) 


def compile_theano_version2(number_of_matrices): 
    assert number_of_matrices > 0 
    L = [tt.matrix() for _ in xrange(number_of_matrices)] 
    concatenated_L = tt.concatenate(L, axis=0) 
    res = tt.dot(concatenated_L.T, concatenated_L) 
    return theano.function(L, res) 


def compile_theano_version3(): 
    concatenated_L = tt.matrix() 
    res = tt.dot(concatenated_L.T, concatenated_L) 
    return theano.function([concatenated_L], res) 


def numpy_version1(*L): 
    assert len(L) > 0 
    n = L[0].shape[1] 
    res = np.zeros((n, n), dtype=L[0].dtype) 
    for M in L: 
     res += np.dot(M.T, M) 
    return res 


def numpy_version2(*L): 
    concatenated_L = np.concatenate(L, axis=0) 
    return np.dot(concatenated_L.T, concatenated_L) 


def main(): 
    iteration_count = 100 
    number_of_matrices = 20 
    n = 300 
    min_x = 400 
    dtype = 'float64' 
    theano_version1 = compile_theano_version1(number_of_matrices, n, dtype) 
    theano_version2 = compile_theano_version2(number_of_matrices) 
    theano_version3 = compile_theano_version3() 
    L = [np.random.standard_normal(size=(x, n)).astype(dtype) 
     for x in range(min_x, number_of_matrices + min_x)] 

    start = timeit.default_timer() 
    numpy_res1 = np.sum(numpy_version1(*L) 
         for _ in xrange(iteration_count)) 
    print 'numpy_version1', timeit.default_timer() - start 

    start = timeit.default_timer() 
    numpy_res2 = np.sum(numpy_version2(*L) 
         for _ in xrange(iteration_count)) 
    print 'numpy_version2', timeit.default_timer() - start 

    start = timeit.default_timer() 
    theano_res1 = np.sum(theano_version1(*L) 
         for _ in xrange(iteration_count)) 
    print 'theano_version1', timeit.default_timer() - start 

    start = timeit.default_timer() 
    theano_res2 = np.sum(theano_version2(*L) 
         for _ in xrange(iteration_count)) 
    print 'theano_version2', timeit.default_timer() - start 

    start = timeit.default_timer() 
    theano_res3 = np.sum(theano_version3(np.concatenate(L, axis=0)) 
         for _ in xrange(iteration_count)) 
    print 'theano_version3', timeit.default_timer() - start 

    assert np.allclose(numpy_res1, numpy_res2) 
    assert np.allclose(numpy_res2, theano_res1) 
    assert np.allclose(theano_res1, theano_res2) 
    assert np.allclose(theano_res2, theano_res3) 


main() 

Quando viene eseguito questo stampe (qualcosa di simile)

numpy_version1 1.47830819649 
numpy_version2 1.77405482179 
theano_version1 1.3603150303 
theano_version2 1.81665318145 
theano_version3 1.86912039489 

Il passaggio afferma, mostrando che la Theano e le versioni numpy entrambi calcolano lo stesso risultato ad alto grado di precisione. Chiaramente questa precisione si riduce se si utilizza float32 anziché float64.

I risultati di temporizzazione mostrano che l'approccio vettoriale non può essere preferibile, dipende dalle dimensioni della matrice. Nell'esempio sopra le matrici sono grandi e l'approccio di non concatenazione è più veloce ma se i parametri n evengono modificati nella funzione main per essere molto più piccoli, l'approccio di concatenazione è più veloce. Altri risultati potrebbero rimanere in esecuzione su una GPU (solo versioni Theano).

Answer 3

Oltre a @ di DanielRenshaw risposta, se aumentiamo il numero di matrici a 1000, la funzione compile_theano_version1 sarà resa RuntimeError: maximum recursion depth exceeded e compile_theano_version2 sembra richiedere un'eternità per la compilazione.

C'è una correzione per questo utilizzando typed_list:

def compile_theano_version4(number_of_matrices, n): 
    import theano.typed_list 
    L = theano.typed_list.TypedListType(tt.TensorType(theano.config.floatX, broadcastable=(None, None)))() 
    res, _ = theano.scan(fn=lambda i: tt.dot(L[i].T, L[i]), 
         sequences=[theano.tensor.arange(number_of_matrices, dtype='int64')]) 
    return theano.function([L], res.sum(axis=0)) 

Inoltre, ho impostato il tipo di dati di tutte le variabili rilevanti per float32 e corsi @ sceneggiatura di DanielRenshaw su GPU, si è scoperto che @ suggerimento di Divakar (theano_version3) è il più efficiente in questo caso. Anche se come ha detto @DanielRenshaw, l'utilizzo di una matrice enorme potrebbe non essere sempre una buona pratica.

Di seguito sono riportate le impostazioni e le uscite sulla mia macchina.

iteration_count = 100 
number_of_matrices = 200 
n = 300 
min_x = 20 
dtype = 'float32' 
theano.config.floatX = dtype 


numpy_version1 5.30542397499 
numpy_version2 3.96656394005 
theano_version1 5.26742005348 
theano_version2 1.76983904839 
theano_version3 1.03577589989 
theano_version4 5.58366179466