如何基于某个等价关系从矩阵列表中删除重复项？

Question

给定一些对称整数矩阵列表，我希望在以下等价关系下删除所有重复项：

如果存在一些置换s（在{1,...,k}上），使得对于{1,...,k}中的所有i和j，我们有M1_ij = M2_s(i)s(j)，即如果我可以通过同时排列其行和列来从一个矩阵获得另一个矩阵，则两个k x k矩阵M1和M2是等价的。

不幸的是，我的朴素方法（在构建列表时，检查新矩阵的任何置换是否已经在列表中）证明速度太慢。

我能想到的一些可能更快的替代方法是将所有矩阵放入列表中，将它们排列成某种“规范置换”，然后按照这里所述的方式删除重复项。但是，我不确定如何在代码中实现这样的“标准置换”。

为了进一步缩小范围：矩阵相对较小（k <= 4），列表将包含大约5或6位数字的矩阵，并且矩阵的dtype必须是某种整数类型（目前为intc，但我可以更改）。

最终列表的顺序无关紧要，每个等价类的代表都可以生存。整个过程可能需要一些小时，但不需要几天。

有没有一种相对高效的方法来实现这一点？我是否（又一次）错过了一些很酷的NumPy或SciPy设施，可以帮助我处理这个问题？

如要求所示，以下是一些小例子，以演示等价关系是如何工作的：

矩阵{{1, 1, 1}, {1, 2, 0}, {1, 0, 3}}和{{1, 1, 1}, {1, 3, 0}, {1, 0, 2}}是等价的，因为置换{1,2,3}->{1,3,2}将一个转换为另一个。

矩阵{{1, 1, 1}, {1, 2, 0}, {1, 0, 3}}和{{1, 1, 0}, {1, 2, 1}, {0, 1, 3}}不等价，您不能更改1的位置而不进行对角线置换。

- Baum mit Augen

---

@jotasi 是的，我们对每个条目的行索引和列索引使用相同的排列；例如，每个对角线条目将再次成为对角线条目。 - Baum mit Augen

@MSeifert 添加了一些示例，它们足够清晰吗？ - Baum mit Augen

只是一个快速的备注：相同类别的两个矩阵必须具有相同的特征值。这可以作为第一次测试。鉴于您的矩阵规模较小，这可能是相当有效的。 - PAb

说到快速测试必要条件，你可以从比较 np.sort(np.diag(M)) 开始。 - jotasi

为什么不趁机将所有条目排序并进行比较 np.all(np.sort(M.ravel()), np.sort(N.ravel()))？！ - Paul Brodersen

Answer 1

这是一个代数解。我怀疑应该有一个更令人满意的组合解。

如果存在置换矩阵P使得M'=P^{-1} M P，则你说两个矩阵M和M'是等价的。

让我们使用M和M'的特征分解：

M = Q^{-1} D Q

M = Q'^{-1} D' Q'

其中D和D'是包含特征值的对角矩阵，Q和Q'是正交矩阵。

我们可以将相等关系重写为：

D = D'（即这两个矩阵应该具有相同的特征值）

Q' = PQ

测试第二个条件很容易。鉴于Q是正交的，它相当于检查矩阵点积(Q, Q'.T)是否是置换矩阵，即它是否在每行和每列上只有一个“1”。

因此，草案算法如下：

- 取M和M' - 计算M和M'的特征分解(Q, D)和(Q', D')（使用np.linalg.eigh） - 如果它们没有相同的特征值（当然要考虑数值精度），则它们不等价 - 否则，计算np.dot(Q, Q'.T)，并测试它是否是置换矩阵

我认为瓶颈是特征分解，但您只需要对每个矩阵执行一次。希望第一个测试能够快速丢弃许多矩阵。

希望这有所帮助。

- PAb

---

谢谢，我会尝试这个和图表的方法，可能结合一些“桶解决方案”来避免O(n^2)。 - Baum mit Augen

即使您有一种解决小于O(n^2)的等价类的方法，但在某些时候，您仍然需要比较两个矩阵，这将始终是O(n^2)（除非您说您不需要100％的确定性并且只比较少量元素）。 - Jürg W. Spaak

Answer 2

您可以将矩阵视为表示图的邻接/权重矩阵，然后测试这两个图是否同构。 networkx有一个方便的函数（可以通过pip安装）。

import numpy as np
import networkx as nx
from networkx.algorithms.isomorphism import numerical_edge_match

# create matrices
n = 4
a = np.random.randint(0, 10, size=(n,n))
a = a + a.T # i.e. symmetric
b = np.rot90(a, k=2) # i.e. a rotated by 180 degrees
c = np.ones((n,n), dtype=np.int) # counter-example

# create graphs
ga = nx.from_numpy_matrix(a)
gb = nx.from_numpy_matrix(b)
gc = nx.from_numpy_matrix(c)

# test if isomorphic
print "a isomorphic with b:", nx.is_isomorphic(ga, gb, edge_match=numerical_edge_match('weight', 1)) # True
print "a isomorphic with c:", nx.is_isomorphic(ga, gc, edge_match=numerical_edge_match('weight', 1)) # False

- Paul Brodersen

---

看起来很有趣，谢谢。让我试一试，看看能否使用它来摆脱我目前的O(n^2)方法；如果不行，我就得测量一下这是否足够快。 - Baum mit Augen

还可以看看 graph-tool (https://graph-tool.skewed.de/) 和 igraph，它们的核心函数是用 C 和 C++ 写的，因此可能更快。 - Paul Brodersen

Answer 3

只需使用您的规范方法。 在矩阵中搜索最大的条目，将其放置在右上角。 然后根据其条目对第一列和第一行进行排序。

A = np.array([[1,2,3,5],
     [3,6,2,6],
     [3,5,7,2],
     [1,3,6,3]])
a = np.where(A == np.amax(A))
sort_colums = np.argsort(A[a[0]].ravel())[::-1]
sort_rows = np.argsort(A[:,a[1]].ravel())[::-1]
Col_sorted = A[:,sort_colums]
Equiv_class = Col_sorted[sort_rows]
#returns [[7, 5, 3, 2],
          [6, 3, 1, 3],
          [3, 2, 1, 5],
          [2, 6, 3, 6]]

正如评论中指出的那样，这仅适用于矩阵条目仅出现一次的情况。如果它们出现多次但不是那么频繁，那么可以通过生成多个等价类来调整此方法。

- Jürg W. Spaak

---

矩阵中不一定需要单个最大元素，我认为这很重要，不是吗？（虽然我还不熟悉那些函数，但我会去查看。） - Baum mit Augen

这确实不适用于像 [[1,1,1],[1,1,0],[1,0,1]] 这样的矩阵。"TypeError: only length-1 arrays can be converted to Python scalars"。 - Baum mit Augen

@BaummitAugen 实际上，它只适用于具有一个特定最大值的矩阵。但是对于这些矩阵，算法的复杂度为O(n log(n))，可以将矩阵转换为等价类。 我更新了代码（我犯了一个小错误） - Jürg W. Spaak