Python NumPy向量化

Question

我正在尝试编写未加权顶点覆盖问题的列表右启发式代码。背景如下：

顶点覆盖问题：在顶点覆盖问题中，我们给出一个无向图G =（V，E），其中V是顶点集合，E是边集合。我们需要找到最小的V'子集，使得V'覆盖G。如果集合V'覆盖图G，则表示集合V'中的所有顶点都至少与图中的一条边相连。

列表右启发式：该算法非常简单。给定一个顶点列表V = [v1，v2，... vn]，其中n是G中的顶点数，如果i>j且vi和vj在图G中由一条边连接，则称vi为vj的右邻居。我们初始化一个覆盖C = {}（空集）并从右到左扫描V。在任何时候，假设当前正在扫描的顶点是u。如果u有至少一个不在C中的右邻居，则将u添加到c中。整个V只被扫描一次。

我正在同时解决多个图形（具有相同的顶点但不同的边缘）的问题。

我用Python编写了List Right Heuristic。我能够将其向量化以同时解决多个图形，但我无法向量化原始for循环。我使用邻接矩阵表示图形。我想知道它是否可以进一步向量化。以下是我的代码：

def list_right_heuristic(population: np.ndarray, adj_matrix: np.ndarray):
    adj_matrices = np.matlib.repmat(adj_matrix,population.shape[0], 1).reshape((population.shape[0], *adj_matrix.shape))

    for i in range(population.shape[0]):
        # Remove covered vertices from the graph. Delete corresponding edges
        adj_matrices[i, np.outer(population[i], population[i]).astype(bool)] = 0

    vertex_covers = np.zeros(shape=population.shape, dtype=population.dtype)
    for index in range(population.shape[-1] - 1, -1, -1):
        # Get num of intersecting elements (for each row) in right neighbors and vertex_covers
        inclusion_rows = np.sum(((1 - vertex_covers) * adj_matrices[..., index])[..., index + 1:], axis=-1).astype(bool)
        # Only add vertices to cover for rows which have at least one right neighbor not in vertex cover
        vertex_covers[inclusion_rows, index] = 1

    return vertex_covers

我有p个图需要同时解决，其中p = population.shape[0]。每个图具有相同的顶点但不同的边缘。population数组是一个二维数组，其中每一行表示已在覆盖中的图G的顶点。我只尝试找到不在覆盖中的顶点。因此，为此原因，将所有行和列设置为覆盖中的顶点为0，即删除相应的边缘。现在，算法应该理论上仅返回不在覆盖中的顶点。 所以在第一个循环中，我只是将邻接矩阵中的相应行和列设置为0（行和列中的所有元素都将为零）。接下来，我从右到左遍历顶点的二维数组，并找出每一行中不在vertex_covers中的右邻居的数量。为此，我首先找到未被覆盖的顶点（1-vertex_covers），然后将其乘以相应的adj_matrices列（或行，因为adj矩阵是对称的），以获取我们正在扫描的顶点的邻居。然后，我加总了它右侧的所有元素。如果该值大于0，则至少有一个不在vertex_covers中的右邻居。 我对一个这样做得对吗？ 是否有任何方法可以向量化第二个循环（或第一个循环），或者总体加速代码？对于具有1000个以上顶点的大型图形，在一些其他代码中调用此函数数千次。任何帮助都将不胜感激。

- Aditya369

---

population和adj_matrix的形状是相同的，对吗？而且它们都将是方形数组吗？ - Divakar

邻接矩阵是一个方阵，但顶点集合不一定是。例如，考虑使用邻接矩阵[[0,0,1,0],[0,0,0,1],[1,0,0,1],[0,1,1,0]]编码的图形中的两个不同的顶点覆盖。假设覆盖1包含第2个和第4个顶点，而覆盖2包含第2个和第3个顶点。那么顶点集合将是[[0,1,0,1], [0,1,1,0]]。虽然通过添加更多带有随机覆盖的行可以轻松地使顶点集合成为方阵。 - Aditya369

Answer 1

您可以使用np.einsum在索引之间执行许多复杂操作。在您的情况下，第一个循环可以这样执行：

adj_matrices[np.einsum('ij, ik->ijk', population, population).astype(bool)] = 0

我花了一些时间才理解einsum的工作原理。我发现这个Stack Overflow问题非常有帮助。

顺便说一下，你的代码给了我以下的语法错误：

SyntaxError: can use starred expression only as assignment target

我不得不重新编写函数的第一行：

adj_matrices = np.matlib.repmat(adj_matrix,population.shape[0], 
        1).reshape((population.shape[0],) + adj_matrix.shape)

- Ramon Crehuet