没有 fill_diagonal 的 inf 的 Python 矩阵对角线

Question

我需要将矩阵的对角线元素设置为 Inf。

一个简单的方法是使用np.fill_diagonal。

np.fill_diagonal(my_matrix, float('inf')

但是，fill_diagonal 会修改输入矩阵，而不是返回填充对角线的新矩阵。 这对我不起作用。我需要填充对角线而不修改原始矩阵。

当然，我可以克隆原始矩阵，所以我将始终保留原始矩阵的副本。但是我不太喜欢这个解决方案，因为我会经常更新原始矩阵，因此每次需要对角线为 inf 时我都必须复制它。

是否有一个函数可以执行与 fill_diagonal 相同的操作，但不修改输入矩阵？像这样的东西:

new_matrix = np.fill_diagonal(original_matrix, float('inf') 

为什么我需要这个:

我的矩阵是点之间的距离矩阵，我想在每一步计算两个最近的点。当然这个矩阵的对角线是0(因为一个点到它自身的距离是0)。因此，为了确保不采取相同的点，我的解决方案是将对角线设置为 Inf。

但是，一旦找到这两个点，我需要计算这两点与其余点之间的距离的平均值，因此我实际上需要对角线为 0 而不是 Inf。

目前我正在做的是:

• 用 Inf 填充对角线
• 找到 2 个最近的点
• 用 0 填充对角线

• 计算这两点与其余点之间的平均距离。

  # fill diagonal with Inf to avoid taking the diagonals
  np.fill_diagonal(data, float('inf'))  
  # find the minimum distance
  idx = np.argmin(data)
  # fill the diagonals back to 0
  np.fill_diagonal(data, 0.0) 
  # get the coordinates of the minimum distance
  row, col =  np.unravel_index(idx,data.shape)
  # compute the new node as the average distance between the two points
  new_node = np.mean((data[:,row],data[:,col]),0)
  # replace the first node (row) with the new node
  data[:,row] = new_node
  data[row,:] = new_node.T
  # delete the second node (col) from the matrix
  data = np.delete(data, col, 0)  # delete row
  data = np.delete(data, col, 1)  # delete column
  

但是我不喜欢将对角线设置为 Inf 然后再设置回 0 的想法，我更愿意将一个函数传递给 argmax ，该函数返回对角线填充 Inf 的数据，而不实际修改矩阵数据。

类似于:

idx = np.argmin(return_filled_diagonals(data, float('Inf'))
# here I can operate with data as usual since it has not been modified.

Answer 1

方法#1

您正在寻找的魔法就在 NumPy strides 中这使我们能够看到没有对角线元素的数组 View ，因此不再占用内存空间。这是获得这样的 View 的实现 -

def nodiag_view(a):
    m = a.shape[0]
    p,q = a.strides
    return np.lib.stride_tricks.as_strided(a[:,1:], (m-1,m), (p+q,q))

让我们看一下示例运行来验证它是否是一个 View -

In [124]: a  # Input array
Out[124]: 
array([[ 0, 61, 43, 26, 21],
       [20,  0, 78, 29, 64],
       [34, 49,  0, 64, 60],
       [36, 96, 67,  0, 75],
       [36, 85, 40, 74,  0]])

# Get the no-diag view
In [125]: a_nodiag = nodiag_view(a)

# Lets's verify by changing elements in the view and that should change
# elements in the original array too
In [126]: a_nodiag[:] = 999

In [127]: a
Out[127]: 
array([[  0, 999, 999, 999, 999],
       [999,   0, 999, 999, 999],
       [999, 999,   0, 999, 999],
       [999, 999, 999,   0, 999],
       [999, 999, 999, 999,   0]])

最后，让我们看看如何设置它来解决您的整个问题 -

def argmin_without_diag(a):
    a_nodiag = nodiag_view(a)
    idx_nodiag = np.argmin(a_nodiag)
    m = a.shape[0]
    idx = idx_nodiag + np.unravel_index(idx_nodiag, (m-1,m))[0]+1
    return  np.unravel_index(idx, a.shape)

示例运行 -

In [142]: a
Out[142]: 
array([[ 0, 60, 79, 55, 77],
       [62,  0, 86, 84, 25],
       [32, 96,  0, 74, 89],
       [24, 33, 64,  0, 93],
       [14, 74, 30, 44,  0]])

In [143]: argmin_without_diag(a)
Out[143]: (4, 0)

<小时/>

方法#2

如果你同时担心内存和性能，可以暂时将对角线设置为infnite，然后获取argmin索引，然后放回原始对角线值。因此，实际上输入数组没有改变。实现看起来像这样 -

def argmin_without_diag_replacement(a):
    # Store diagonal values
    vals = a.ravel()[::a.shape[1]+1].copy()

    # Set diag ones as infinites
    a.ravel()[::a.shape[1]+1] = np.inf

    # Get argmin index
    idx = np.argmin(a)

    # Put back the original diag values
    a.ravel()[::a.shape[1]+1] = vals
    return np.unravel_index(idx, a.shape)

因此，对于 (n x n) 形状的数组，临时数组将只有 n 个元素。

示例运行 -

In [237]: a
Out[237]: 
array([[  0.,  95.,  57.,  75.,  92.],
       [ 37.,   0.,  69.,  71.,  62.],
       [ 42.,  72.,   0.,  30.,  57.],
       [ 41.,  80.,  94.,   0.,  26.],
       [ 36.,  45.,  71.,  76.,   0.]])

In [238]: argmin_without_diag_replacement(a)
Out[238]: (3, 4)

运行时测试

In [271]: a = np.random.randint(11,99,(1000,1000)).astype(float)

In [272]: np.fill_diagonal(a,0)

In [273]: %timeit argmin_without_diag(a)
1000 loops, best of 3: 1.76 ms per loop

In [274]: %timeit argmin_without_diag_replacement(a)
1000 loops, best of 3: 688 µs per loop