从 scipy CSR 矩阵索引到 numpy 数组的最有效方法？

Posted 2023-03-12

【中文标题】从 scipy CSR 矩阵索引到 numpy 数组的最有效方法？

【英文标题】：Most efficient way to index into a numpy array from a scipy CSR matrix?

【发布时间】：2018-07-01 14:19:35

【问题描述】：

我有一个形状为(4000, 3) 的numpy ndarray X，其中X 中的每个样本都是一个3D 坐标（x，y，z）。

我有一个形状为(4000, 4000) 的scipy csr 矩阵nn_rad_csr，它是从sklearn.neighbors.radius_neighbors_graph(X, 0.01, include_self=True) 生成的最近邻图。

nn_rad_csr.toarray()[i] 是一个形状 (4000,) 稀疏向量，其二进制权重（0 或 1）与节点 X[i] 的最近邻图中的边相关联。

例如，如果nn_rad_csr.toarray()[i][j] == 1 则X[j] 在X[i] 的最近邻居半径内，而0 的值表示它不是邻居。

我想做的是有一个函数radius_graph_conv(X, rad)，它返回一个数组Y，它是X，由其邻居的值平均。我不确定如何利用 CSR 矩阵的稀疏性来有效地执行radius_graph_conv。我在下面有两个简单的 graph conv 实现。

import numpy as np
from sklearn.neighbors import radius_neighbors_graph, KDTree

def radius_graph_conv(X, rad):
    nn_rad_csr = radius_neighbors_graph(X, rad, include_self=True)
    csr_indices = nn_rad_csr.indices
    csr_indptr  = nn_rad_csr.indptr
    Y = np.copy(X)
    for i in range(X.shape[0]):
        j, k = csr_indptr[i], csr_indptr[i+1]
        neighbor_idx = csr_indices[j:k]
        rad_neighborhood = X[neighbor_idx] # ndim always 2
        Y[i] = np.mean(rad_neighborhood, axis=0)
    return Y

def radius_graph_conv_matmul(X, rad):
    nn_rad_arr = radius_neighbors_graph(X, rad, include_self=True).toarray()
    # np.sum(nn_rad_arr, axis=-1) is basically a count of neighbors

    return np.matmul(nn_rad_arr / np.sum(nn_rad_arr, axis=-1), X)

有没有更好的方法来做到这一点？使用 knn 图，它是一个非常简单的功能，因为邻居的数量是固定的，您可以只索引 X，但是对于基于半径或密度的最近邻居图，您必须使用 CSR，（或如果您使用的是 kd 树，则为数组）。

【问题讨论】：

您获取每行非零索引的方法是较快的方法之一。另一种选择是转换为lil 格式。 rows 属性的每个元素都是同一个行列表。

另一种方法是使用X[csr_indices] 一次索引所有邻居，将cumsum 应用于此，然后使用csr_indptr 计算子和均值。

【参考方案1】：

这是利用 csr 格式的直接方法。您的 matmul 解决方案可能在幕后做了类似的事情。但是我们通过利用它是一个邻接矩阵来保存一个查找（来自.data 属性）；此外，diffing .indptr 应该比对等量的求和更有效。

>>> import numpy as np
>>> from scipy import sparse
>>> 
# create mock data
>>> A = np.random.random((100, 100)) < 0.1
>>> A = (A | A.T).view(np.uint8)
>>> AS = sparse.csr_matrix(A)
>>> X = np.random.random((100, 3))
>>> 
# dense solution for reference
>>> Xa = A @ X / A.sum(axis=-1, keepdims=True)
# sparse solution
>>> XaS = np.add.reduceat(X[AS.indices], AS.indptr[:-1], axis=0) / np.diff(AS.indptr)[:, None]
>>> 
# check they are the same
>>> np.allclose(Xa, XaS)
True

【讨论】：

这很棒。不幸的是，我使用的框架没有像np.add.reduceat 这样的高级索引减少函数，而且我必须使用支持的函数（matmul）进行梯度链接，所以我选择了密集的解决方案。感谢您提供有关区分 indptr 的提示！