使用 3D numpy 数组的直观方式

Question

让我们举一个简单的例子，我有数据数组

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

现在假设我有一组数据，在多个时间步长中采用二维数组的形式。为简单起见，再次假设此数据只是 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 数组来使其工作，但现在第二个维度充当 "time" 索引，这有点违反直觉，因为在AA 它是表示时间的第一个维度...理想情况下，我想要一个解决方案，其中 AA 和 BB 都在第三个维度中具有时间索引！

编辑：

因为 dot(a, b)[i,j,k,m] = sum(a[i,j,:] * b[k,:,m]) 来自文档，我想知道我想要实现的目标是否不可能？这似乎很奇怪，因为这应该是人们可能想要的一种相对普遍的东西......

Answer 1

我不确定我是否同意您对 "intuitive" 的定义 - 对我来说，用数组的第一维表示时间索引似乎更自然。由于 numpy 数组对于每个 2x2 子矩阵是 row-major by default, this ordering of the dimensions will give you better 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 用作 shorthand for np.matmul:

来使其更加紧凑

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

Answer 2

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]: Q@A@Q.T
Out[102]: 
array([[ 1.75,  2.75],
       [ 4.  ,  4.5 ]])
In [103]: Q@AA@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)

使用 3D numpy 数组的直观方式

Intuitive way to use 3D numpy arrays

python

numpy

matrix-multiplication