Зацикливание (или векторизация) матриц переменной длины в Theano

Question

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

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

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

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

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

Редактировать:

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

Редактировать:

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

Результаты синхронизации показывают, что векторизованный подход не может быть предпочтительным, он зависит от размеров матрицы. В приведенном выше примере матрицы большие, и подход без конкатенации быстрее, но если n а также 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) является наиболее эффективным в этом случае. Хотя, как сказал @DanielRenshaw, использование огромной матрицы не всегда может быть хорошей практикой.

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

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