Python NumPy 向量化

Question

我正在尝试为未加权的顶点覆盖问题编写所谓的 List Right Heuristic 代码。背景如下:

顶点覆盖问题:在顶点覆盖问题中，我们得到一个无向图 G = (V, E)，其中 V 是顶点集，E 是边集。我们需要找到最小的集合 V'，它是 V 的子集，使得 V' 覆盖 G。如果图中的所有边在 V' 中至少有一个顶点，则称集合 V' 覆盖图 G。/p>

List Right Heuristic:算法非常简单。给定顶点列表 V = [v1, v2, ... vn] 其中 n 是 G 中的顶点数，如果 i > j 并且 vi 和 vj 由边连接，则称 vi 是 vj 的右邻居在图 G 中。我们启动一个覆盖 C = {}(空集)并从右到左扫描 V。在任何时候，假设正在扫描的当前顶点是 u。如果你至少有一个右邻居不在 C 中，那么你就被添加到 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 中已经在封面中的顶点。我只是想找到不在封面中的顶点。因此，出于这个原因，将 cover 中所有顶点的行和列设置为 0，即，我正在删除相应的边。启发式算法理论上应该现在只返回不在封面中的顶点。 所以在第一个 for 循环中，我只是将邻接矩阵中相应的行和列设置为 0(行和列中的所有元素都将为零)。接下来，我将从右到左遍历二维顶点数组，并在不在 vertex_covers 中的每一行中找到右邻居的数量。为此，我首先找到不在 cover (1 - vertex_covers) 中的顶点，然后将其与 adj_matrices 中的相应列相乘(或行，因为 adj 矩阵是对称的)以获得我们正在扫描的那个顶点的邻居。然后我将所有元素加到这个右边。如果该值大于 0，则至少有一个右邻不在 vertex_covers 中。 我这样做正确吗？ 有什么方法可以向量化第二个 for 循环(或第一个 for 循环)或总体上加速代码？在一些其他代码中为大型图(具有 1000 多个顶点)调用此函数数千次。任何帮助，将不胜感激。

Answer 1

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

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

我花了一些时间才明白 einsum 是如何工作的。我找到了 this SO question很有帮助。

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

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)