如何使用 NumPy 中的列表列表矢量化高级索引？

Question

以下代码在使用纯 Python 时运行时间为 45 秒。

for iteration in range(maxiter):
    for node in range(n):
        for dest in adjacency_list[node]:
            rs[iteration + 1][dest] += beta * rs[iteration][node] / len(adjacency_list[node])

但是，通过简单地将 rs 初始化为 numpy ndarray 而不是 python 列表，代码在 145 秒内运行。我真的不知道为什么 numpy 在这个数组索引上花费 3 倍的时间。

我的想法是尽可能多地矢量化，但只设法矢量化了 beta/len(adjacency_list[node]) 的乘法。此代码在 77 秒内运行。

beta_over_out_degree = np.array([beta / len(al) for al in adjacency_list])
for iteration in range(1, maxiter + 1):
    r_next = np.full(shape=n, fill_value=(1 - beta) / n)
    f = beta_over_out_degree * r
    for i in range(n):
        r_next[adjacency_list[i]] += f[i]

    r = np.copy(r_next)
    rs[iteration] = np.copy(r)

问题在于adjacency_list 是具有不同列大小的列表的列表，具有 100 000 行和 1-15 列。 具有邻接矩阵的更标准方法（至少作为正常的 ndarray）不是一种选择，因为对于 n=100 000，它的 (n,n) 形状太大而无法分配给内存。

有没有办法使用它的索引来向量化 numpy 高级索引（也许把它变成一个 numpy ndarray）？

我也非常感谢任何其他速度提示。 提前致谢！

编辑：感谢@stevemo，我设法创建了具有csr_matrix 功能的adjacency_matrix，并将其用于迭代乘法。程序现在只需 2 秒即可运行！

for iteration in range(1, 101):
    rs[iteration] += rs[iteration - 1] * adjacency_matrix

Answer 1

如果我理解正确，这可以通过使用邻接矩阵的矩阵幂的单线公式来完成。

根据您的原始代码 sn-p，您似乎有一些 n 节点网络，邻接信息存储为 adjacency 中的列表列表，并且您有一个与每个关联的值 r节点在迭代k+1 时的值是beta 乘以它在迭代k 的每个邻居的r 之和。 （您的循环以相反的方向构造它，但同样的事情。）

如果你不介意将你的adjacency list-of-lists 改造成更标准的adjacency matrix，这样A_ij = 1 如果ij 是邻居，则为0，那么你可以完成内部两个循环使用简单的矩阵乘积，r[k+1] = beta * (A @ r[k])。

并遵循这个逻辑，r[k+2] = beta * (A @ (beta * (A @ r[k]))) = (beta * A)**2 @ r[k] 或一般来说，

r[k] = (beta * A)**k @ r[0]

让我们在一个小型网络上试试这个：

# adjacency matrix
A = np.array([
    [0, 1, 1, 0, 0],
    [1, 0, 1, 0, 0],
    [1, 1, 0, 1, 0],
    [0, 0, 1, 0, 1],
    [0, 0, 0, 1, 0]
])

# initial values
n = 5
beta = 0.5
r0 = np.ones(n)
maxiter = 10

# after one iteration
print(beta * (A @ r0))
# [1.  1.  1.5 1.  0.5]

# after 10 iterations
print(np.linalg.matrix_power((beta * A), maxiter) @ r0)
# [2.88574219 2.88574219 3.4921875  1.99414062 0.89257812]