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