python - Python 中的快速、小型和重复矩阵乘法

Question

我正在寻找一种使用 Python/Cython/Numpy 快速将许多 4x4 矩阵相乘的方法，有人可以提供任何建议吗？

为了展示我目前的尝试，我有一个需要计算的算法

A_1 * A_2 * A_3 * ... * A_N 

每一个地方

A_i != A_j

它在 Python 中的一个例子:

means = array([0.0, 0.0, 34.28, 0.0, 0.0, 3.4])
stds = array([ 4.839339, 4.839339, 4.092728, 0.141421, 0.141421, 0.141421])

def fn():
    steps = means+stds*numpy.random.normal(size=(60,6))
    A = identity(4)
    for step in steps:
        A = dot(A, transform_step_to_4by4(step))
%timeit fn()

1000 loops, best of 3: 570 us per loop

在 Cython/Numpy 中实现此算法比使用 Eigen/C++ 进行所有优化的等效代码慢大约 100 倍。不过，我真的不想使用 C++。

Answer 1

如果您必须调用 Python 函数来生成每个要相乘的矩阵，那么您的性能基本上就完蛋了。但是，如果您可以向量化 transform_step_to_4by4 函数，并让它返回一个形状为 (n, 4, 4) 的数组，那么您可以使用 matrix_multiply< 节省一些时间:

import numpy as np
from numpy.core.umath_tests import matrix_multiply

matrices = np.random.rand(64, 4, 4) - 0.5

def mat_loop_reduce(m):
    ret = m[0]
    for x in m[1:]:
        ret = np.dot(ret, x)
    return ret

def mat_reduce(m):
    while len(m) % 2 == 0:
        m = matrix_multiply(m[::2], m[1::2])
    return mat_loop_reduce(m)

In [2]: %timeit mat_reduce(matrices)
1000 loops, best of 3: 287 us per loop

In [3]: %timeit mat_loop_reduce(matrices)
1000 loops, best of 3: 721 us per loop

In [4]: np.allclose(mat_loop_reduce(matrices), mat_reduce(matrices))
Out[4]: True

您现在有 log(n) 个 Python 调用而不是 n，有利于 2.5 倍的加速，对于 n = 1024 将接近 10 倍。显然 matrix_multiply 是一个 ufunc，并且因此有一个 .reduce 方法，这将允许您的代码在 Python 中运行无循环。我一直无法让它运行，不断收到神秘错误:

In [7]: matrix_multiply.reduce(matrices)
------------------------------------------------------------
Traceback (most recent call last):
  File "<ipython console>", line 1, in <module>
RuntimeError: Reduction not defined on ufunc with signature