量子で解くはずが、最小費用流で解けてしまった ── 治験マッチングを JijModeling 2 と OMMX で書いて踏んだ6つの罠

量子アニーリングで解くつもりで定式化した問題が、最小費用流で厳密に解けてしまいました。多項式時間です。量子を持ち出す理由が、最初の版の時点で消えていました。

題材は治験の患者マッチングです。NEDO の量子懸賞金事業で、複数の治験に患者をまとめて割り当てる問題を QUBO にしています。モデルは JijModeling 2 で書き、OMMX を経由して SCIP・OpenJij などのソルバに同じインスタンスを渡します。一つのモデルから複数のソルバに出せるので、「古典で解けるのでは？」という問いに同じ土俵で答えられる。それが狙いでした。

その過程で六つの罠を踏みました。どれもエラーにならず、もっともらしい数字が出る種類のものです。順に並べます。

モデル

変数は「適格性を満たす (患者, 治験) の組」ごとに一つです。患者は高々一つの治験に入り、治験には定員があります。目的は、候補ごとのスコアの和から、同じ治験に同じ層（年齢層など）の患者が偏った分を引いたものです。

import jijmodeling as jm

@jm.Problem.define("TrialMatching", sense=jm.ProblemSense.MAXIMIZE)
def trial_matching(problem: jm.DecoratedProblem):
    w = problem.Float("w", ndim=1)
    n = w.len_at(0, latex="n")
    lam = problem.Float("lam", latex=r"\lambda")
    x = problem.BinaryVar("x", shape=(n,))

    PT = problem.Placeholder("PT", ndim=1,
        dtype=(jm.DataType.NATURAL, jm.DataType.NATURAL))  # (患者, 候補)
    TT = problem.Placeholder("TT", ndim=1,
        dtype=(jm.DataType.NATURAL, jm.DataType.NATURAL))  # (治験, 候補)
    SS = problem.Placeholder("SS", ndim=1,
        dtype=(jm.DataType.NATURAL, jm.DataType.NATURAL))  # 同一治験・同一層の候補対
    cap = problem.Natural("cap", ndim=1)
    n_p = problem.Natural("n_p")
    n_t = cap.len_at(0, latex="n_t")

    problem += jm.sum(w[i] * x[i] for i in n)
    problem -= lam * jm.sum(x[e[0]] * x[e[1]] for e in SS)

    problem += problem.Constraint("patient",
        (jm.sum(x[e[1]] for e in PT if e[0] == p) <= 1
         for p in jm.range(n_p)))
    problem += problem.Constraint("capacity",
        (jm.sum(x[e[1]] for e in TT if e[0] == t) <= cap[t]
         for t in jm.range(n_t)))

疎な構造なので、患者と候補の対応は (患者番号, 候補番号) のタプルの列で持ち、ストリームを if で絞って総和を取っています。この書き方が罠6につながります。

罠1：定員制約をどう書いても、問題は簡単なまま

最初の版には、二次項がありませんでした。層の偏りは候補ごとの重みに織り込めばよいと考えていたからです。

ところが「患者は高々1つ」「治験は高々 k」「目的は線形」という組合せは、二部 b-マッチングそのものです。最小費用流に帰着して、多項式時間で厳密に解けます。患者12人・治験2本・8シードで測ると、最小費用流は全シードで厳密解と一致しました。貪欲法でも一致しました。

定員  λ     厳密    QUBO    最小費用流  貪欲   流れ損失
  2  0.0   16.77   16.77   16.77    16.77   0.00
  5  0.0   37.30   37.30   37.30    37.30   0.00

難しさを生むのは、同じ治験に入る患者同士の関係を見る二次項です。層が揃っているかどうかは、個々の患者の属性では決まりません。だから本質的に二次で、最小費用流には載りません。

定員  λ     厳密    QUBO    最小費用流  貪欲   流れ損失
  2  1.5   15.51   15.51   14.52    14.52   0.99
  3  1.5   19.64   19.64   18.05    18.05   1.59
  5  1.5   22.18   22.18   19.11    19.11   3.07

損失は定員とともに開きます（0.99 → 1.59 → 3.07）。定員が大きいほど、同じ治験の中で層が衝突する組が増えるからです。逆に言うと、定員1〜3の小さい設定で測ると二次項の効果を見誤ります。

教訓：定式化を変えるたびに、最小費用流（か、その問題に合う多項式アルゴリズム）を一緒に回す。QUBO が厳密解を返して喜ぶ前に、古典の多項式解法が同じ値を返していないかを確かめる。

罠2：「ちょうど k 人」のペナルティが、他の制約を食い破る

JijModeling に移す前、QUBO を手で組んでいた版では、定員を「ちょうど k 人」の二乗ペナルティで書いていました。

P_cap * (Σ_{i∈T} x_i − k)^2

これは二つの理由で壊れました。

一つ目は単純で、候補が定員に満たないシードでは「ちょうど k」を満たせません。6シード中4つで実行可能解が存在しませんでした。

二つ目のほうが厄介です。展開すると、対角に P_cap × (1 − 2k) が出ます。k = 3 なら −5 × P_cap です。この大きな負の対角項が、「同じ患者を二つの治験に入れない」ためのペナルティ P_conflict = max(w) + 1 を上回ってしまいました。あるシードでは、QUBO の最良解の目的値が 35.06 で、「高々 k」で求めた本当の最適値は 30.80 でした。制約を破った解が、最小エネルギーになっていたのです。

「高々 k」を二進のスラック変数で表す形に変えて、8シードすべてで厳密解と一致しました。JijModeling と OMMX に移してからは、このスラック変数を自分で書く必要がなくなっています。

教訓：二乗ペナルティは展開してから係数を見る。対角項は k に比例して大きくなり、別の制約のペナルティより大きくなることがある。

罠3：厳密解より良い値が出たら、まず疑う

アニーリング系のソルバは、制約をペナルティとして目的関数に足して解きます。ペナルティ係数は小さいほうが解の質が良くなりがちです。そのぶん、制約を破った解のほうがエネルギーが低くなることもあります。

返ってきたサンプルの目的値をそのまま比べると、厳密な最適値を上回る数字が出ます。最大化問題で厳密解より良い値が出るはずはありません。出たなら、それは定員を超えて詰め込んだ解か、一人の患者を二つの治験に入れた解です。

これを二度踏みました。一度目に気づいて直し、ソルバを足したときにもう一度やりました。

for smp in samples:
    sel = selected(smp)
    if not is_feasible(cands, trials, sel):   # 先に絞る
        continue
    best = max(best, objective(sel))

教訓：目的値を計算する前に、必ず実行可能性で絞る。厳密解より良い値は、成果ではなく不具合の兆候。

罠4：OpenJij は、Q に出てこない変数を結果から落とす

自分で組んだ QUBO を OpenJij に直接渡した版で踏みました。結果の record.sample を列の位置で読んでいたら、選ばれた候補がずれました。

原因は、Q に一度も現れない変数は、結果の変数一覧から落ちることです。係数がすべて 0 になったスラックのビットなどが該当します。一つ落ちると、それより後ろの列が全部一つずつずれます。エラーは出ません。

# 列は位置ではなく変数名で引く
col = {v: k for k, v in enumerate(res.variables)}
for smp in res.record.sample:
    sel = [i for i in range(n_x) if i in col and smp[col[i]]]

教訓：サンプルの列は res.variables を通して変数名で引く。変数の個数が入力と同じだと決めつけない。

罠5：OMMX 経由の SCIP では、変数名が「’0′, ‘1’, …」になっている

SCIP が探索中に見つけた解を全部集めようとして、getSols() の各解から x_0、x_1…という名前で値を引きました。結果は、どの解も何も選んでいないように見えました。

OMMX のアダプタが SCIP に渡す変数の名前は、OMMX の変数ID を文字列にしたもの（'0'、'1'、…）でした。JijModeling で付けた x という名前ではありません。存在しない名前で引いたので、空の選択が静かに返っていました。

m = OMMXPySCIPOptAdapter(instance).solver_input
m.optimize()
varmap = {v.name: v for v in m.getVars()}
for sol in m.getSols():
    sel = [i for i in range(n_x)
           if str(i) in varmap and m.getSolVal(sol, varmap[str(i)]) > 0.5]

しかも「選ばれた候補がゼロ」は、制約を全部満たす実行可能解です。罠3の実行可能性チェックも、そのまますり抜けます。

教訓：ソルバ側の変数名は、アダプタを通したあとで一度 getVars() を出力して確かめる。モデルで付けた名前が残っているとは限らない。

罠6：絞り込んだストリームで、制約検出が型エラーになる（2.8.0）

JijModeling 2.8.0 では、上のモデルを eval() すると E-TE0017 で落ちました。if e[0] == p で絞ったストリームの要素の型が ElementOf[...] のまま簡約されず、e[1] の添字付けに失敗していました。引っかかっていたのは制約検出の処理でした。

回避策は、制約検出を切ることです。

instance = trial_matching.eval(data, constraint_detection=False)  # 2.8.0 の回避策

ただし、これには副作用の心配がありました。制約検出は、SOS1 のような構造を見つけて古典ソルバに渡す機能です。切ったまま比べると、SCIP を不当に弱く見せているおそれがあります。量子と古典を比べる実験で、それは致命的です。

JIJ の Discord（#jijmodeling_日本語）に報告したところ、その日のうちに再現・原因特定・修正まで進み、2日後に 2.9.0 として出ました。2.9.0 で制約検出を有効に戻して測り直し、検出の有無で SCIP の結果が変わらないことを確かめています。

教訓：回避策で機能を切ったら、その機能が比較の公平性に効いていないかを、直ったあとに測り直して確かめる。

まとめ

六つとも、例外は投げません。罠1は「量子が厳密解を出した」、罠2と罠3は「厳密解より良い値が出た」、罠4と罠5は「解がずれた／空だった」、罠6は「回避策で動いた」。どれも、それらしい結果として通り過ぎていくものでした。

JijModeling と OMMX の組み合わせは、こうした確認を回すのに向いています。一つのモデルから SCIP にも OpenJij にも出せるので、「古典ではどうか」を毎回同じインスタンスで測れます。罠1に気づけたのも、最小費用流と QUBO と厳密解を、同じ表に並べていたからでした。