使用3D numpy的陣列

直觀的方式讓作爲簡單的例子，其中我有數據陣列使用3D numpy的陣列

A = np.asarray([[1,3], [2,4]])

而這個數據是以下一個簡單的變換將被變換爲另一種形式：

Q = np.asarray([[-0.5,1], [1,0.5]]) 
B = np.dot(Q,np.dot(A,Q.T)) 
print B

現在假設我有一組數據，它採用2d數組的形式幾個時間步長。爲了簡單起見，再次假定這個數據只是A被複制3個時間步驟。我們可以將此數據表示爲尺寸爲(2,2,N)的3d陣列，其中N =3在這種情況下。第三維然後代表數據的時間索引。現在很自然的要求有一種簡單的方式來轉換上面的數據，但是對於每個時間步，通過直觀的3D數組相乘，但是我只能做出以下非直觀的工作：

# Create the 3d data array 
AA = np.tile(A,(3,1,1)) # shape (3,2,2) 
BB = np.dot(Q,np.dot(AA,Q.T)) 

print np.all(BB[:,0,:] == B) # Returns true

所以用這個方法我沒有重鑄Q陣列，使其工作，但是現在的第二個維度充當這是一個有點反直覺的，因爲在AA這是第一個「時間」指數表示時間的維度...理想情況下，我想要一個解決方案，其中AA和BB在第三維中具有時間索引！

編輯：

由於dot(a, b)[i,j,k,m] = sum(a[i,j,:] * b[k,:,m])從文檔我想知道如果我想實現是不可能的？這似乎很奇怪，因爲這應該是一個相對普遍的事情，人們可能會希望...

來源

2016-09-26 Jack

隨着'numpy'的'第一軸是最外面的一個，所以'（N，2,2）'形狀可能更加自然，例如由np.array（[A1，A2，A3]）'產生，即3個2x2陣列的連接。但是這種'點'旋轉也可以與（2，N，2）或（2,2，N）一起使用。 'dot'通常適用於最後一次暗淡，第二次暗淡。 'np.matmul'和'np.einsum'是替代品。 – hpaulj

@hpaulj那麼如何在這種情況下執行數組乘法？ – Jack

In [91]: A=np.array([[1,3],[2,4]]) 
In [92]: Q=np.array([[-.5,1],[1,.5]]) 
In [93]: B=np.dot(Q,np.dot(A,Q.T)) 
In [94]: B 
Out[94]: 
array([[ 1.75, 2.75], 
     [ 4. , 4.5 ]])

相同的計算與einsum：

In [95]: np.einsum('ij,jk,kl',Q,A,Q) 
Out[95]: 
array([[ 1.75, 2.75], 
     [ 4. , 4.5 ]])

如果我做的A多個副本 - 一個新的第一維：

In [96]: AA = np.array([A,A,A]) 
In [97]: AA.shape 
Out[97]: (3, 2, 2) 
... 
In [99]: BB=np.einsum('ij,pjk,kl->pil',Q,AA,Q) 
In [100]: BB 
Out[100]: 
array([[[ 1.75, 2.75], 
     [ 4. , 4.5 ]], 

     [[ 1.75, 2.75], 
     [ 4. , 4.5 ]], 

     [[ 1.75, 2.75], 
     [ 4. , 4.5 ]]])

BB有（3,2,2）的形狀。

的新望matmul（@運算符）讓我做同樣的事情

In [102]: [email protected]@Q.T 
Out[102]: 
array([[ 1.75, 2.75], 
     [ 4. , 4.5 ]]) 
In [103]: [email protected]@Q.T 
Out[103]: 
array([[[ 1.75, 2.75], 
     [ 4. , 4.5 ]], 

     [[ 1.75, 2.75], 
     [ 4. , 4.5 ]], 

     [[ 1.75, 2.75], 
     [ 4. , 4.5 ]]])

隨着einsum它也很容易與最後一維的工作：

In [104]: AA3=np.stack([A,A,A],-1) # newish np.stack 
In [105]: AA3.shape 
Out[105]: (2, 2, 3) 
In [106]: np.einsum('ij,jkp,kl->ilp',Q,AA3,Q) 
Out[106]: 
array([[[ 1.75, 1.75, 1.75], 
     [ 2.75, 2.75, 2.75]], 

     [[ 4. , 4. , 4. ], 
     [ 4.5 , 4.5 , 4.5 ]]]) 
In [107]: _.shape 
Out[107]: (2, 2, 3)

來源

2016-09-26 19:28:48 hpaulj

非常好 - 我只是在試用'einsum'！看起來這是一個非常靈活的總結方式。你有沒有關於你的頭腦的知識如何比較它的性能和點？ – Jack

另外還有一件事 - 你的第一個解決方案稍微不正確，它應該是'BB = np.einsum'（'ij，pjk，lk-> pil'，Q，AA，Q）'儘管在這種情況下它不會標記問題因爲'Q'是對稱的！ – Jack

確實;我在心裏做了一個註釋，說'Q.T'需要切換索引順序，但沒有驗證它。 – hpaulj

我不確定我是否同意你的「直覺」的定義 - 對我來說，第一個代表時間索引似乎更自然數組的維數。由於numpy數組默認情況下爲row-major，因此所有維度的排序將爲每個2x2子矩陣提供更好的locality of reference，因爲所有元素都將位於相鄰的存儲器地址中。

儘管如此，有可能通過使用np.matmul和調換每個中間陣列，以適應您的例子來爲(2, 2, N)陣列工作：

CC = np.repeat(A[..., None], 3, -1) # shape (2, 2, 3) 
DD = np.matmul(Q.T, np.matmul(CC.T, Q)).T 

print(DD.shape) 
# (2, 2, 3) 

print(repr(DD)) 
# array([[[ 1.75, 1.75, 1.75], 
#   [ 2.75, 2.75, 2.75]], 

#  [[ 4. , 4. , 4. ], 
#   [ 4.5 , 4.5 , 4.5 ]]])

在Python 3.5+可以讓這個更緊湊的使用@ operator作爲簡寫np.matmul：

DD = (Q.T @ (CC.T @ Q)).T

來源

2016-09-26 19:13:17

使用3D numpy的陣列

回答

相關問題