在Theano中循环(或矢量化)可变长度矩阵

Question

我有一个矩阵列表L,其中每个项目M都是一个x*n矩阵(x是一个变量,n是一个常量).

我想计算以下Python代码所做的M'*M所有项目的总和L(M'是转置M):

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

实际上我想在Theano中实现它(它不支持可变长度的多维数组),我不想将所有矩阵填充到相同的大小,因为这会浪费太多空间(一些矩阵可能非常大) ).

有一个更好的方法吗？

编辑:

L 在Theano汇编之前就已经知道了.

编辑:

收到了来自@DanielRenshaw和@Divakar的两个优秀答案,情绪上难以选择接受.

Answer 1

鉴于在Theano编译需要发生之前已知矩阵的数量,人们可以简单地使用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和min_x参数更改main为更小,则串联方法更快.在GPU上运行时可能会保留其他结果(仅限Theano版本).

Answer 2

您可以沿着第一个轴填充输入数组，即所有x' 的总和。因此，我们最终会得到一个高(X,n)数组，其中X =x1+x2+x3+..... 这可以被转置，并且它与其自身的点积将是 shape 的所需输出(n,n)。所有这一切都是通过利用强大的点积的纯矢量化解决方案来实现的。因此，实施将是 -

\n\n

# Concatenate along axis=0\nLcat = np.concatenate(L,axis=0)\n\n# Perform dot product of the transposed version with self\nout = Lcat.T.dot(Lcat)\n

\n\n

运行时测试并验证输出 -

\n\n

In [116]: def vectoized_approach(L):\n     ...:   Lcat = np.concatenate(L,axis=0)\n     ...:   return Lcat.T.dot(Lcat)\n     ...: \n     ...: def original_app(L):\n     ...:   n = L[0].shape[1]\n     ...:   res = np.zeros((n,n))\n     ...:   for M in L:\n     ...:     res += np.dot(M.T, M)\n     ...:   return res\n     ...: \n\nIn [117]: # Input\n     ...: L = [np.random.rand(np.random.randint(1,9),5)for iter in range(1000)]\n\nIn [118]: np.allclose(vectoized_approach(L),original_app(L))\nOut[118]: True\n\nIn [119]: %timeit original_app(L)\n100 loops, best of 3: 3.84 ms per loop\n\nIn [120]: %timeit vectoized_approach(L)\n1000 loops, best of 3: 632 \xc2\xb5s per loop\n

\n