如何使用 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，我设法创建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倍数它的邻居在 iter k。（你的循环在相反的方向构造它，但同样的事情。）

如果您不介意将您的adjacency列表列表改造成更标准的邻接矩阵，例如A_ij = 1ifij是邻居，否则为 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]