numpy中的大点矩阵数组乘法

Question

给定两个大的numpy数组，一个用于3D点列表，另一个用于转换矩阵列表。假设两个列表之间存在1对1的对应关系，我正在寻找最佳方法来计算由相应矩阵变换的每个点的结果数组。

我的解决方法做，这是（在下面的示例代码中看到“TEST4”），它工作得很好用小数组使用切片，但是失败，因为我的方法是如何内存浪费是:)

大型阵列

import numpy as np 
COUNT = 100 
matrix = np.random.random_sample((3,3,)) # A single matrix 
matrices = np.random.random_sample((COUNT,3,3,)) # Many matrices 
point = np.random.random_sample((3,))  # A single point 
points = np.random.random_sample((COUNT,3,)) # Many points 

# Test 1, result of a single point multiplied by a single matrix 
# This is as easy as it gets 
test1 = np.dot(point,matrix) 
print 'done' 

# Test 2, result of a single point multiplied by many matrices 
# This works well and returns a transformed point for each matrix 
test2 = np.dot(point,matrices) 
print 'done' 

# Test 3, result of many points multiplied by a single matrix 
# This works also just fine 
test3 = np.dot(points,matrix) 
print 'done' 

# Test 4, this is the case i'm trying to solve. Assuming there's a 1-1 
# correspondence between the point and matrix arrays, the result i want 
# is an array of points, where each point has been transformed by it's 
# corresponding matrix 
test4 = np.zeros((COUNT,3)) 
for i in xrange(COUNT): 
    test4[i] = np.dot(points[i],matrices[i]) 
print 'done' 

用一个小阵列，这工作正常。对于大数组，（COUNT = 1000000）测试＃4有效，但速度相当慢。

有没有办法让Test＃4更快？假定不使用循环？

Answer 1

您可以使用numpy.einsum。这里有5点矩阵和5点的例子：

In [49]: matrices.shape 
Out[49]: (5, 3, 3) 

In [50]: points.shape 
Out[50]: (5, 3) 

In [51]: p = np.einsum('ijk,ik->ij', matrices, points) 

In [52]: p[0] 
Out[52]: array([ 1.16532051, 0.95155227, 1.5130032 ]) 

In [53]: matrices[0].dot(points[0]) 
Out[53]: array([ 1.16532051, 0.95155227, 1.5130032 ]) 

In [54]: p[1] 
Out[54]: array([ 0.79929572, 0.32048587, 0.81462493]) 

In [55]: matrices[1].dot(points[1]) 
Out[55]: array([ 0.79929572, 0.32048587, 0.81462493]) 

以上是做matrix[i] * points[i]（即乘以右侧），但我只是重读的问题，发现你的代码使用points[i] * matrix[i]。你可以做到这一点通过切换einsum的指标和参数：

In [76]: lp = np.einsum('ij,ijk->ik', points, matrices) 

In [77]: lp[0] 
Out[77]: array([ 1.39510822, 1.12011057, 1.05704609]) 

In [78]: points[0].dot(matrices[0]) 
Out[78]: array([ 1.39510822, 1.12011057, 1.05704609]) 

In [79]: lp[1] 
Out[79]: array([ 0.49750324, 0.70664634, 0.7142573 ]) 

In [80]: points[1].dot(matrices[1]) 
Out[80]: array([ 0.49750324, 0.70664634, 0.7142573 ]) 

Answer 2

拥有多个变换矩阵没有多大意义。您可以组合变换矩阵作为this question：

如果我想申请矩阵A，然后是B，然后是C，我要加倍矩阵以相反的顺序np.dot(C,np.dot(B,A))

所以，你可以通过预先计算是节省一些内存空间矩阵。然后将一组矢量应用于一个变换矩阵应该很容易处理（理由）。

我不知道为什么你需要在一百万个载体上进行一百万次转换，但我建议购买一个更大的RAM。

编辑： 没有办法减少操作，没有。除非你的变换矩阵具有稀疏性，对角性等特定属性，否则你将不得不运行所有乘法和求和。但是，处理这些操作的方式可以跨核心进行优化和/或在GPU上使用矢量操作。

另外，python速度很慢。您可以尝试使用NumExpr跨核心分割numpy。或者在C++上使用BLAS框架（特别是快速;））