在 Theano 中循环（或向量化）可变长度矩阵

我有一个矩阵列表L，其中每个项目M is a x*n矩阵 （x是一个变量，n是一个常数）。

我想计算总和M'*M对于所有项目L (M'是转置M）如以下 Python 代码所示：

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

实际上我想在Theano中实现这个（它不支持可变长度多维数组），并且我不想将所有矩阵填充到相同的大小，因为这会浪费太多空间（有些矩阵可能非常大） ）。

有一个更好的方法吗？

Edit:

L在 Theano 编译之前就已知。

Edit:

从@DanielRenshaw 和@Divakar 收到了两个出色的答案，在情感上很难选择一个接受。

鉴于在需要进行 Theano 编译之前矩阵的数量是已知的，因此可以简单地使用常规的 Python 矩阵列表。

这是一个完整的示例，显示了 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 and min_x参数改变在main函数要小得多，然后串联方法更快。在 GPU 上运行时可能会出现其他结果（仅限 Theano 版本）。