在 np.matrix 中使用索引条件进行操作的最佳方法答案

【问题标题】：Best way to do operations in np.matrix with conditions on the indexes在 np.matrix 中使用索引条件进行操作的最佳方法
【发布时间】：2022-11-01 11:29:35
【问题描述】：

我正在寻找在取决于索引条件的 numpy 矩阵中执行操作的最佳方法。

我正在处理的矩阵是一个对称方阵，特别是它是一个加权邻接矩阵。

目前我有三个嵌套循环，计算成本很高。

下面的代码记录了我如何执行操作以及循环执行期间的条件。

# matrix is a numpy.matrix square matrix, in particular a weighted adjacency matrix
result_vector = []
for i in range(matrix.shape[0]):
    aux = 0
    for j in range(matrix.shape[0]):
         if j != i:
            for k in range(matrix.shape[0]):
                if k != j:
                    aux += (matrix[i,j]*matrix[i,k])*(1 - matrix[j,k])
    result_vector.append(aux)
result_vector = np.array(result_vector)

我尝试使用numpy.einsum，但由于操作中的减法我没有成功。

有没有办法在避免循环的同时执行操作？

【问题讨论】：

首先，最好坚持使用普通的 numpy 数组；不鼓励使用np.matrix，因为它通常比有用更令人困惑。 @ 是矩阵乘法运算符。替换您的if 逻辑可能很棘手，至少在没有彻底可视化正在发生的事情的情况下并非如此。制作一个或多个 mask 数组可能有助于在对角线上/在对角线之外为真或假。替换循环需要考虑对整个数组的操作，而不是逐个元素。

标签： python numpy loops numpy-ndarray

【解决方案1】：

您与einsum 走在了正确的轨道上。您可以在一次调用中将该函数与您想要的任何数组一起使用，这意味着减法可以临时存储为新数组并在einsum 中使用。您需要小心以某种方式考虑einsum 遍历所有索引的事实，因此必须手动强制执行您的约束。

def func2_einsum(matrix):
    # Create helper array
    one_minus_matrix = 1-matrix

    # For the subexpression matrix[i,k])*(1 - matrix[j,k])
    # Save the result for any i, j, k
    T = np.einsum('ik,jk->ij', matrix, one_minus_matrix)
    # Save the result for j == k
    T_j_eq_k = np.einsum('ij,jj->ij', matrix, one_minus_matrix)
    # Enforce j != k
    T = T - T_j_eq_k
    # Enforce i != j
    np.fill_diagonal(T, 0)

    # Compute & return the full expression matrix[i,j] * T[i,j]
    return np.einsum('ij,ij->i', matrix, T)

对基于循环的方法进行的一些基准测试表明，einsum 方法的速度提高了 3 个数量级（取决于 N，数组的大小）。

N = 3
A = np.random.rand(N, N)

%timeit func1_loops(A)
15.6 us ± 248 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
%timeit func2_einsum(A)
12.1 us ± 64.9 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

N = 10
A = np.random.rand(N, N)

%timeit func1_loops(A)
708 us ± 8.41 us per loop (mean ± std. dev. of 7 runs, 1000 loops each)
%timeit func2_einsum(A)
13.2 us ± 172 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

N = 100
A = np.random.rand(N, N)

%timeit func1_loops(A)
792 ms ± 2.64 ms per loop (mean ± std. dev. of 7 runs, 1 loops each)
%timeit func2_einsum(A)
439 us ± 1.99 ns per loop (mean ± std. dev. of 7 runs, 1000 loops each)

【讨论】：