Python SciPy optimize.quadratic_assignment用法及代码示例

本文简要介绍 python 语言中 scipy.optimize.quadratic_assignment 的用法。

用法:

scipy.optimize.quadratic_assignment(A, B, method='faq', options=None)#

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

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

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

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

图匹配试图最大化相同的目标函数。该算法可以被认为是找到两个图的节点的对齐方式，以最小化引起的边不一致的数量，或者，在加权图的情况下，是平方边权重差的总和。

请注意，二次分配问题是NP-hard。这里给出的结果是近似值，不能保证是最优的。

参数 ：：

A： 二维阵列，正方形

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

B： 二维阵列，正方形

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

method： {‘faq’, ‘2opt’} 中的字符串(默认值：‘faq’)

用于解决问题的算法。 ‘faq’(默认)和‘2opt’可用。

options： 字典，可选

求解器选项字典。所有求解器都支持以下内容：

最大化 布尔(默认值：假)

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

partial_match 二维整数数组，可选(默认值：无)

修复部分匹配。也称为“seed” [2]。

每一行partial_match指定一对匹配的节点：nodepartial_match[i, 0]的A与节点匹配partial_match[i, 1]的B.数组有形状(m, 2)，其中m不大于节点数，\(n\).

rng {无，整数， numpy.random.Generator ，

numpy.random.RandomState}, optional

如果种子是无(或np.random)， 这numpy.random.RandomState使用单例。如果种子是一个 int，一个新的RandomState使用实例，播种种子.如果种子已经是一个Generator或者RandomState实例然后使用该实例。

对于 method-specific 选项，请参阅 show_options('quadratic_assignment') 。

返回 ：：

res： OptimizeResult

OptimizeResult 包含以下字段。

col_ind 一维数组

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

乐趣 浮点数

解决方案的客观值。

尼特 int

优化期间执行的迭代次数。

注意：

默认方法‘faq’使用Fast Approximate QAP算法[1]；它通常提供速度和准确性的最佳组合。方法 ‘2opt’ 的计算成本可能很高，但可能是一种有用的替代方法，或者它可以用于细化其他方法返回的解决方案。

参考：

[1]

J.T. Vogelstein, J.M. Conroy, V. Lyzinski, L.J. Podrazik, S.G. Kratzer, E.T.哈雷，D.E.菲什金德，R.J. Vogelstein 和 C.E. Priebe，“用于图形匹配的快速近似二次规划”，PLOS one，vol。 10，没有。第 4 页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”、模式识别。 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_ind 和 fun 之间的关系，使用 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