SCIPに整数性のあるLPを解かせてみる。公式のPySCIPOptを使う。
from itertools import combinations
from pyscipopt import Model, quicksum, SCIP_PARAMSETTING
n = 4
model = Model("graph_matching")
xs = {
(i, j): model.addVar(lb=0.0, name=f"x_{i}_{j}")
for i, j in combinations(range(n), 2)
}
for v in range(n):
model.addCons(quicksum(x for (i, j), x in xs.items() if i == v or j == v) == 1)
model.setObjective(quicksum(xs.values()), "minimize")
model.optimize()
sol = {(i, j): model.getVal(x) for (i, j), x in xs.items()}
print(f"objective: {model.getObjVal()}")
print(f"solution: {sol}")
制約は$n=4$の完全グラフに関する完全マッチング多面体を表していて、変数が6つ、等号制約が4つなので、基底解はゼロになる変数を少なくとも2つ含む。例えば$(1, 0, 0, 0, 0, 1)$など。しかし、SCIPは基底解を返さない:
objective: 2.0
solution: {(0, 1): 0.3333333333333333, (0, 2): 0.3333333333333333, (0, 3): 0.3333333333333333, (1, 2): 0.33333333333333337, (1, 3): 0.33333333333333337, (2, 3): 0.33333333333333337}
SCIPはLPソルバではないので、分枝限定法の1回目のLPを行う前に問題を変形したりカットを追加したりはしているはずで、この結果はおかしくはない。ただ、READMEを眺めてpresolveをオフにしてみたり、他の設定を弄ったりしても特に結果に変わりなかった:
model.setPresolve(SCIP_PARAMSETTING.OFF)
model.setHeuristics(SCIP_PARAMSETTING.OFF)
model.disablePropagation()
model.writeProblem("graph_matching.lp")でLPを書き出して直接SCIPに解かせてみることにする。LPファイルの中身は以下のようになる:
\ SCIP STATISTICS
\ Problem name : graph_matching
\ Variables : 6 (0 binary, 0 integer, 0 implicit integer, 6 continuous)
\ Constraints : 4
Minimize
Obj: +1 x01 +1 x02 +1 x03 +1 x12 +1 x13 +1 x23
Subject to
c1: +1 x01 +1 x02 +1 x03 = +1
c2: +1 x01 +1 x12 +1 x13 = +1
c3: +1 x02 +1 x12 +1 x23 = +1
c4: +1 x03 +1 x13 +1 x23 = +1
Bounds
End
これをSCIPのシェルを直接使って解いてもやはり結果は変わらない:
SCIP> read graph_matching.lp
SCIP> set presolving emphasis off
SCIP> set heuristics emphasis off
SCIP> set separating emphasis off
SCIP> optimize
SCIP> display solution
objective value: 2
x01 0.333333333333333 (obj:1)
x02 0.333333333333333 (obj:1)
x03 0.333333333333333 (obj:1)
x12 0.333333333333333 (obj:1)
x13 0.333333333333333 (obj:1)
x23 0.333333333333333 (obj:1)
write transproblemでpresolve後のLPを出力してみる:
\ SCIP STATISTICS
\ Problem name : t_graph_matching.lp
\ Variables : 6 (0 binary, 0 integer, 0 implicit integer, 6 continuous)
\ Constraints : 12
Minimize
Obj: +1 t_x01 +1 t_x02 +1 t_x03 +1 t_x12 +1 t_x13 +1 t_x23
Subject to
c1: +1 t_x01 +1 t_x02 +1 t_x03 = +1
c2: +1 t_x01 +1 t_x12 +1 t_x13 = +1
c3: +1 t_x02 +1 t_x12 +1 t_x23 = +1
c4: +1 t_x03 +1 t_x13 +1 t_x23 = +1
SSTcut_0_1: -1 t_x01 +1 t_x02 <= +0
SSTcut_0_2: -1 t_x01 +1 t_x03 <= +0
SSTcut_0_3: -1 t_x01 +1 t_x12 <= +0
SSTcut_0_4: -1 t_x01 +1 t_x13 <= +0
SSTcut_0_5: -1 t_x01 +1 t_x23 <= +0
SSTcut_1_2: -1 t_x02 +1 t_x03 <= +0
SSTcut_1_3: -1 t_x02 +1 t_x12 <= +0
SSTcut_1_4: -1 t_x02 +1 t_x13 <= +0
Bounds
0 <= t_x01 <= 1
0 <= t_x02 <= 1
0 <= t_x03 <= 1
0 <= t_x12 <= 1
0 <= t_x13 <= 1
0 <= t_x23 <= 1
End
SSTcutが犯人で、$x_{01} \ge x_{02}$のような制約を等号にすると確かに上のような「基底解」が出ることになる。[1] 問題はset presolve emphasis offでもこのカットが無効にならないことだが、とりあえずset misc usesymmetry 0とすることで解決した:
SCIP> display solution
objective value: 2
x03 1 (obj:1)
x12 1 (obj:1)
後から気づいたが、PySCIPOptからも同等の設定がmodel.setIntParam("misc/usesymmetry", 0)でできるようだ。
結論として、set presolve emphasis offやそれに対応するPySCIPOptのsetPresolve(SCIP_PARAMSETTING.OFF)はpresolveのフェーズで起こることを完全に無しにしてくれるわけではない。おそらくmisc/usesymmetry = 0を加えても完全ではなく、stackoverflowでも別の人が別の問題で困っているようだった。
そもそもLPソルバが使いたいならSCIPの背後にあるSoPlexを直接使うべきなんだろうけど、SoPlexには公式のPythonバインディングがないという問題がある。[2] Pythonで単体法のLPソルバが使いたい時はSCIPではなく、CBC/CLPやHiGHSを検討したほうがいいかもしれない。特にCBC/CLPのPythonバインディングであるCyLPは、READMEを見た感じだとLPソルバを使って自前でILPを解くような状況を意識しているようだ。（試してない。）
たぶんこのカットは問題の対称性がある時に等価な枝を調べるのを避けるためにある ↩︎
