如何使用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列表，代码可以在145s内运行。我真的不知道为什么numpy花费3倍的时间进行此数组索引。

我的想法是对尽可能多的向量进行矢量化处理，但是仅设法对beta/len(adjacency_list[node])的乘法进行矢量化处理。该代码在77s内运行。

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是具有不同列大小的列表的列表，具有10万行和1-15列。至少不使用普通的ndarray作为邻接矩阵的更标准方法是不可行的，因为对于n = 100 000，其形状（n，n）太大，无法分配给内存。

是否有任何方法可以使用其索引进行numpy高级索引（可能将其转换为numpy ndarray？]

我也非常感谢其他速度提示。预先感谢！

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

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

Answer 1

如果我对您的理解正确，可以使用邻接矩阵的矩阵幂使用一线公式来完成。

根据您的原始代码段，似乎您有一些n节点网络，其邻接信息存储为adjacency中的列表列表，并且您具有与每个节点相关联的值r，因此迭代k+1的值是beta乘以迭代器r的每个邻居的k之和。 （您的循环以相反的方向构造了这个，但是是相同的。）

[如果您不介意将adjacency列表重组为更标准的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]