在给定向量 a 的情况下，使用 .sum 在矩阵 F(a[i], a[j]) 上执行三角和的有效方法

Question

我有一个向量a，需要对两个索引进行求和，比如

for i in (range, n): 
    for j in (i+1, n):
        F(a[i] - a[j])

其中F 是一个函数：sum 提醒我们对数组的上三角形求和。

我阅读了Fastest way in numpy to sum over upper triangular elements with the least memory 上的有趣帖子并进行了试验：ARRAY.sum 确实是一种对上三角矩阵元素求和的非常快速的方法。

要将方法应用于我的案例，我首先需要定义一个数组，例如

A[i,j] = F(a[i],a[j])

然后计算

(A.sum() - np.diag(A).sum())/2

我当然可以通过两个 for 循环定义数组 A，但我想知道是否有更快、更简单的方法。

在另一种情况下，函数F 简单地等于

F = a[i]*a[j]

我可以写

def sum_upper_triangular(vector):
    A = np.tensordot(vector,vector,0)
    return (A.sum() - np.diag(A).sum())/2

这比直接用 sum() 求和或嵌套 for 循环要快得多。

如果F更清晰，例如

np.exp(a[i] - a[j])

我想知道最有效的方法。

非常感谢

Answer 1

如果我理解正确，您希望执行以下操作：

result = []
n = len(a)
for i in range(n-1): 
    for j in range(i+1, n):
        result.append(F(a[i] - a[j]))

对于某些功能F。此外，矩阵元素之间的运算可以是任何其他运算（例如乘法*）。下面是一种不使用 for 循环的方法：

iu = np.triu_indices(n, k=1)
A_j, A_i = np.meshgrid(a, a)
res = F(A_i[iu] - A_j[iu])  # e.g. F = np.exp

解释（n=5）：

A_i = [[a[0], a[0], a[0], a[0], a[0],
       [a[1], a[1], a[1], a[1], a[1],
       [a[2], a[2], a[2], a[2], a[2],
       [a[3], a[3], a[3], a[3], a[3],
       [a[4], a[4], a[4], a[4], a[4]]

A_j = [[a[0], a[1], a[2], a[3], a[4],
       [a[0], a[1], a[2], a[3], a[4],
       [a[0], a[1], a[2], a[3], a[4],
       [a[0], a[1], a[2], a[3], a[4],
       [a[0], a[1], a[2], a[3], a[4]]

A_i[iu] = [a[0], a[0], a[0], a[0], a[1], a[1], a[1], a[2], a[2], a[3]]
A_j[iu] = [a[1], a[2], a[3], a[4], a[2], a[3], a[4], a[3], a[4], a[4]]

然后执行元素计算并应用元素 F：

F(A_i[iu] - A_j[iu]) = [ 
    F(a[0] - a[1]), F(a[0] - a[2]), F(a[0] - a[3]), F(a[0] - a[4]),
    F(a[1] - a[2]), F(a[1] - a[3]), F(a[1] - a[4]),
    F(a[2] - a[3]), F(a[2] - a[4]),
    F(a[3] - a[4])]

Answer 2

您可以使用scipy.spatial.pdist 并将您想要的任何功能设置为指标。作为奖励，pdist 仅计算非对角三角形，因此您无需将其从 sum 中删除

from scipy.spatial.distance import pdist

def sum_upper_tri(arr, F = lambda x, y: x*y):
    return pdist(arr.reshape(arr.shape[0], -1), metric = F).sum()/2

如果你想要超快的东西，你需要numba:

from numba import jit

@jit
def sum_upper_tri_jit(arr, F = lambda x, y: x * y):
    out = 0
    for i in range(1, len(arr)):
        for j in range(i + 1, len(arr)):
            out += F(arr[i], arr[j])
    return out / 2

还没有找到@njit 的方法，但如果可以的话，它会快得多。

在任何情况下，每个预期 F 的专门构建的函数都会快得多。比如exp(|x-y|)的情况（提醒exp(x-y)不是对称的：x-y != y-x）

from numba import njit

@njit
def sum_upper_tri_exp(arr):
    out = 0
    for i in range(1, len(arr)):
        for j in range(i + 1, len(arr)):
            out += np.exp(np.abs(arr[i] - arr[j]))
    return out / 2

这比上面的快了大约 100 倍

如果你不想求和，你可以使用：

from numba import njit

@njit
def sum_upper_tri_exp(arr):
    out = []
    for i in range(1, len(arr)):
        for j in range(i + 1, len(arr)):
            out += [np.exp(arr[i] - arr[j])]
    return out