如何在 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

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

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

如果您不介意将adjacency列表列表改革为更标准的邻接矩阵，这样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]