python - 在 Python 中对矩阵的选定元素求和

Question

我有一个 [n x n] 矩阵，其中包含属于不同组的值，以及一个 [1 x n] 向量，定义每个元素属于哪个组。 (n通常~1E4，在本例中n=4)

我想计算一个将属于同一组的所有元素相加得到的矩阵。

我使用 np.where() 来计算每个组的元素所在的索引。 当我使用计算出的索引时，我没有获得预期的元素，因为我选择的是位置对而不是范围(我习惯了 Matlab，我可以简单地选择 M(idx1,idx2) )。

import numpy as np

n=4
M = np.random.rand(n,n)
print(M)

# This vector defines to which group each element belong
belongToGroup = np.array([0, 1, 0, 2])

nGroups=np.max(belongToGroup);

# Calculate a matrix obtained by summing elements belonging to the same group
M_sum = np.zeros((nGroups+1,nGroups+1))
for g1 in range(nGroups+1):
    idxG1 = np.where(belongToGroup==g1)
    for g2 in range(nGroups+1):
        idxG2 = np.where(belongToGroup==g2)
        print('g1 = ' + str(g1))
        print('g2 = ' + str(g2))
        print(idxG1[0])
        print(idxG2[0])
        print(M[idxG1[0],idxG2[0]])
        print(np.sum(M[idxG1[0],idxG2[0]]))
        M_sum[g1,g2]=np.sum(M[idxG1[0],idxG2[0]])

print('')
print('Example of the problem:')
print('Elements I would like to sum to obtain M_sum[0,0]')
print(M[0:2,0:2])
print('Elements that are summed instead')
print(M[[0,1],[0,1]])

问题示例: 在上面的示例中，元素 M_sum[0,0] 应该是 M[0,0]、M[0,1]、M[1,0] 和 M[1,1] 的和 相反，它被计算为 M[0,0] 和 M[1,1] 之和

Answer 1

在 MATLAB 中，使用 2 个列表(实际上是矩阵)进行索引会选择一个 block 。另一方面，numpy 尝试相互广播索引数组，并返回选定的点。其行为类似于 sub2ind 在 MATLAB 中的行为。

In [971]: arr = np.arange(16).reshape(4,4)                                      
In [972]: arr                                                                   
Out[972]: 
array([[ 0,  1,  2,  3],
       [ 4,  5,  6,  7],
       [ 8,  9, 10, 11],
       [12, 13, 14, 15]])
In [973]: i1, i2 = np.array([0,2,3]), np.array([1,2,0])                         

使用 2 个相同大小的一维数组进行索引:

In [974]: arr[i1,i2]
Out[974]: array([ 1, 10, 12])

这实际上返回[arr[0,1], arr[2,2], arr[3,0]]，每个匹配索引点一个元素。

但是，如果我将一个索引转换为“列向量”，它会从行中进行选择，而 i2 会从列中进行选择。

In [975]: arr[i1[:,None], i2]                                                   
Out[975]: 
array([[ 1,  2,  0],
       [ 9, 10,  8],
       [13, 14, 12]])

MATLAB 使 block 索引变得容易，而单独访问则更困难。在 numpy 中， block 访问有点困难，尽管底层机制是相同的。

在您的示例中，i1[0] 和 i2[0] 可以是如下数组:

array([0, 2]), array([3])
(2,) (1,)

形状 (1,) 数组也可以与 (2,) 或 (2,1) 数组一起广播。如果 is[0] 是 np.array([0,1,2])(一个无法与以下对象配对的 (3,) 数组)，您的代码将会失败(2,) 数组。但对于 (2,1)，它会生成 (2,3) block 。