如何使用 NumPy 中的列表列表对高级索引进行矢量化?
How to vectorize advanced indexing with list of lists in NumPy?
以下代码在使用纯 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
如果我没理解错的话,这可以通过使用邻接矩阵的矩阵幂的单线公式来完成。
根据您的原始代码片段,您似乎有一些 n
节点网络,邻接信息存储为 adjacency
中的列表列表,并且您有一个值 r
与每个节点相关联,因此它在迭代 k+1
时的值是 beta
乘以它在迭代 k
时每个邻居的 r
之和。 (你的循环以相反的方向构建它,但同样的事情。)
如果您不介意将您的 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]
以下代码在使用纯 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
如果我没理解错的话,这可以通过使用邻接矩阵的矩阵幂的单线公式来完成。
根据您的原始代码片段,您似乎有一些 n
节点网络,邻接信息存储为 adjacency
中的列表列表,并且您有一个值 r
与每个节点相关联,因此它在迭代 k+1
时的值是 beta
乘以它在迭代 k
时每个邻居的 r
之和。 (你的循环以相反的方向构建它,但同样的事情。)
如果您不介意将您的 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]