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框架（特別是快速;））