петля над (или векторизовать) матрица переменной длиной в Теано

Question

У меня есть список матриц L, где каждый элемент M является x*n матрицей (x является переменной, n является константой).

Я хочу, чтобы вычислить сумму M'*M для всех элементов в L (M' транспонированная M) в качестве следующего кода Python делает:

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

На самом деле я хочу, чтобы реализовать это в Теано (который Безразлично 't поддерживает многомерные массивы переменной длины), и я не хочу, чтобы все матрицы были одинаковыми, потому что это будет слишком много пространства (некоторые из матриц могут быть очень большими).

Есть ли лучший способ сделать это?

Edit:

L известно до компиляции Theano.

Edit:

получил две отличные ответы от @DanielRenshaw и @Divakar, эмоционально трудно выбрать один, чтобы принять.

Answer 1

Вы можете просто проложить входные массивы вдоль первой оси, которая является суммой x. Таким образом, мы получим высокий массив (X,n), где X =x1+x2+x3+..... Это может быть транспонировано, и его точечный продукт с его самостью будет желательным выходом формы (n,n). Все это достигается с помощью чистого векторизованного решения, использующего мощный точечный продукт.Таким образом, реализация будет -

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

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

время выполнения тестов и проверить выход -

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

Учитывая, что число матриц известно до того, как должна произойти компиляция Theano, можно просто использовать обычные списки на языке Python из матриц Theano.

Вот полный пример, показывающий разницу между версиями numpy и Theano.

Этот код был обновлен, чтобы включить сравнение с векторизованным подходом @ Divakar, который работает лучше. Для Теано возможны два векторизованных подхода, где Теано выполняет конкатенацию, а один, где numpy выполняет конкатенацию, результат которой затем передается в 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() 

При запуске Печатает (что-то вроде)

numpy_version1 1.47830819649 
numpy_version2 1.77405482179 
theano_version1 1.3603150303 
theano_version2 1.81665318145 
theano_version3 1.86912039489 

утверждает проход, показывая, что Theano и Numpy версии оба вычислить тот же результат высокой степенью точности. Очевидно, что эта точность будет уменьшаться при использовании float32 вместо float64.

Результаты синхронизации показывают, что векторизованный подход может быть не предпочтительным, он зависит от размеров матрицы. В приведенном выше примере матрицы велики, а подход без конкатенации выполняется быстрее, но если параметры и min_x изменены в функции main намного меньше, то подход конкатенации выполняется быстрее. Другие результаты могут выполняться при работе на GPU (только версии Theano).

Answer 3

В дополнение к @ DanielRenshaw Ответим, если увеличить число матриц до 1000, функция compile_theano_version1 даст RuntimeError: maximum recursion depth exceeded , и compile_theano_version2, похоже, навсегда собирает компиляцию.

Там это исправить это с помощью 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)) 

Кроме того, я установить тип данных всех соответствующих переменных float32 и побежал @ сценарий DanielRenshaw на GPU, оказалось, что @ предложение Divakar (в theano_version3) является наиболее эффективным в этом случае. Хотя, как сказал ДаниэльРеншоу, использование огромной матрицы не всегда может быть хорошей практикой.

Ниже приведены настройки и выходы на моей машине.

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