scipy.optimize.

二次分配#

scipy.optimize.quadratic_assignment(A, B, method='faq', options=None)[源代码][源代码]#

近似求解二次分配问题和图匹配问题。

二次分配解决以下形式的问题:

\[\min_P & \ {\ \text{trace}(A^T P B P^T)}\ \mbox{s.t. } & {P \ \epsilon \ \mathcal{P}}\]

其中 \(\mathcal{P}\) 是所有置换矩阵的集合,而 \(A\)\(B\) 是方阵。

图匹配试图 最大化 相同的目標函數。這個算法可以被認為是找到兩個圖的節點對齊,這最小化了誘導邊的不一致數量,或者在加權圖的情況下,最小化了邊權重差的平方和。

注意,二次分配问题是NP难的。这里给出的结果是近似值,不保证是最优的。

参数:
A2-D 数组, 方阵

上述目标函数中的方阵 \(A\)

B2-D 数组, 方阵

目标函数中的方阵 \(B\)

方法str in {‘faq’, ‘2opt’} (默认: ‘faq’)

用于解决问题的算法。可用的有 ‘faq’ <optimize.qap-faq>`(默认)和 :ref:’2opt’ <optimize.qap-2opt>`。

选项dict, 可选

求解器选项的字典。所有求解器都支持以下选项:

最大化bool (默认: False)

如果为 True,则最大化目标函数。

部分匹配整数的二维数组,可选(默认:无)

修复了部分匹配问题。也称为“种子” [2]

partial_match 的每一行指定了一对匹配的节点:A 的节点 partial_match[i, 0]B 的节点 partial_match[i, 1] 匹配。该数组的形状为 (m, 2),其中 m 不大于节点数 \(n\)

rng{None, int,}

如果 seed 是 None(或 np.random),则使用 numpy.random.RandomState 单例。如果 seed 是 int,则使用新的 RandomState 实例,并以 seed 为种子。如果 seed 已经是 GeneratorRandomState 实例,则使用该实例。

对于特定方法的选项,请参见 show_options('quadratic_assignment')

返回:
res优化结果

OptimizeResult 包含以下字段。

col_ind一维数组

对应于找到的节点 B 的最佳排列的列索引。

乐趣浮动

解决方案的目标值。

nit整数

优化过程中执行的迭代次数。

注释

默认方法 ‘faq’ 使用快速近似QAP算法 [1];它通常提供速度和准确性的最佳组合。方法 ‘2opt’ 可能计算成本较高,但可能是一个有用的替代方案,或者可以用于优化另一种方法返回的解决方案。

参考文献

[1]

J.T. Vogelstein, J.M. Conroy, V. Lyzinski, L.J. Podrazik, S.G. Kratzer, E.T. Harley, D.E. Fishkind, R.J. Vogelstein, 和 C.E. Priebe, “快速近似二次规划用于图匹配,” PLOS one, 第10卷, 第4期, p. e0121002, 2015, DOI:10.1371/journal.pone.0121002

[2]

D. Fishkind, S. Adali, H. Patsolic, L. Meng, D. Singh, V. Lyzinski, C. Priebe, “Seeded graph matching”, Pattern Recognit. 87 (2019): 203-215, DOI:10.1016/j.patcog.2018.09.014

[3]

“2-opt,” 维基百科。 https://en.wikipedia.org/wiki/2-opt

示例

>>> import numpy as np
>>> from scipy.optimize import quadratic_assignment
>>> A = np.array([[0, 80, 150, 170], [80, 0, 130, 100],
...               [150, 130, 0, 120], [170, 100, 120, 0]])
>>> B = np.array([[0, 5, 2, 7], [0, 0, 3, 8],
...               [0, 0, 0, 3], [0, 0, 0, 0]])
>>> res = quadratic_assignment(A, B)
>>> print(res)
     fun: 3260
 col_ind: [0 3 2 1]
     nit: 9

要查看返回的 col_indfun 之间的关系,使用 col_ind 来形成找到的最佳排列矩阵,然后评估目标函数 \(f(P) = trace(A^T P B P^T )\)

>>> perm = res['col_ind']
>>> P = np.eye(len(A), dtype=int)[perm]
>>> fun = np.trace(A.T @ P @ B @ P.T)
>>> print(fun)
3260

或者,为了避免显式构建置换矩阵,可以直接置换距离矩阵的行和列。

>>> fun = np.trace(A.T @ B[perm][:, perm])
>>> print(fun)
3260

虽然一般来说不能保证,但 quadratic_assignment 恰好找到了全局最优解。

>>> from itertools import permutations
>>> perm_opt, fun_opt = None, np.inf
>>> for perm in permutations([0, 1, 2, 3]):
...     perm = np.array(perm)
...     fun = np.trace(A.T @ B[perm][:, perm])
...     if fun < fun_opt:
...         fun_opt, perm_opt = fun, perm
>>> print(np.array_equal(perm_opt, res['col_ind']))
True

以下是一个示例,其中默认方法 ‘faq’ 无法找到全局最优解。

>>> A = np.array([[0, 5, 8, 6], [5, 0, 5, 1],
...               [8, 5, 0, 2], [6, 1, 2, 0]])
>>> B = np.array([[0, 1, 8, 4], [1, 0, 5, 2],
...               [8, 5, 0, 5], [4, 2, 5, 0]])
>>> res = quadratic_assignment(A, B)
>>> print(res)
     fun: 178
 col_ind: [1 0 3 2]
     nit: 13

如果准确性很重要,可以考虑使用 ‘2opt’ 来优化解决方案。

>>> guess = np.array([np.arange(len(A)), res.col_ind]).T
>>> res = quadratic_assignment(A, B, method="2opt",
...                            options = {'partial_guess': guess})
>>> print(res)
     fun: 176
 col_ind: [1 2 3 0]
     nit: 17