2025年1月25日土曜日

Pythonでソルバから双対変数を取り出す方法とその落とし穴

以下の最適化問題が与えられたとする。

$$ \begin{array}{lll} \text{maximize} & f(x)\\ \text{subject to} & g_i(x) \leq 0 & \text{for } i = 1, ..., k\\ & x \in \mathbb R^n\\ \end{array} $$

このうちのひとつの制約$g_j(x)\le 0$を選んで、他の制約を固定しつつこの制約の右辺の定数を微小に増減させた時に目的関数値がどう変わるか知りたいことがある。これを制約$g_j(x) \le 0$のシャドウプライスと呼ぶ。つまり、

$$ \begin{array}{llll} L_j(b) :=& \text{maximize} & f(x)\\ & \text{subject to} & g_i(x) \leq 0 & \text{for } i = 1, ..., j-1, j+1, ..., k\\ && g_j(x) \le b \\ && x \in \mathbb R^n\\ \end{array} $$

としたときの$L_j$の$b=0$における傾きがシャドウプライスである。[1] 微分可能な場合は微分だし、そうでない場合も劣微分を考えればよい。ちなみにLPの場合には$L_j$は区分線形凹関数になることが知られている。[2] これは不等式制約の例だが、等式制約についても同様に定義できる。

LPや一部の凸最適化問題など、ラグランジュ双対問題との強双対性が成り立つものについては、シャドウプライスは制約に対応する双対変数の値に一致する。[3] たいていのソルバは求解の過程で双対問題も同時に解いているので問題が解ければこの値も簡単に得られる。例えばSoPlexをCLIで起動する場合はsoplex <input file> -yなどで参照できる。

Pythonのソルバインターフェースやモデリングツールも普通は双対変数を取り出す機能を備えている。例えばPuLPだったら制約オブジェクトのpiというアトリビュートにシャドウプライスが格納される。[4]

import pulp as pl

problem = pl.LpProblem("problem", pl.LpMaximize)
x = pl.LpVariable("x", lowBound=0)
y = pl.LpVariable("y", lowBound=0)
constr = x + 2 * y <= 1
problem += constr
problem += 2 * x + y
problem.solve()
print(constr.pi)

上の例を実行すると、$x + 2y \le 1$のシャドウプライスが2.0であると正しく表示される。

しかし、「符号の正しいシャドウプライス」を得るためには注意が必要である。一般にシャドウプライスは制約の右辺の増分に対して定義される。これは$x+2y \le 1$ と$1 \ge x + 2y$ではシャドウプライスの符号が逆になると解釈できる。もちろん、シャドウプライスは普通は$f(x) \le 0$とか$Ax \ge b$のように制約のある種の「標準形」に対して定義されるので、非標準的な立式に対して符号がどちらになるべきかは曖昧と言えるかもしれない。ただ、少なくとも使っているツールのシャドウプライスの符号に一貫的でわかりやすいルールがあることは期待したくなる。しかし、現実にはそれはかなわないようだ。

import pulp as pl

x = pl.LpVariable("x", lowBound=0)
y = pl.LpVariable("y", lowBound=0)


def solve(constr: pl.LpConstraint) -> float:
    problem = pl.LpProblem("problem", pl.LpMaximize)
    problem += constr
    problem += 2 * x + y

    problem.solve()
    return constr.pi


print(solve(x + 2 * y <= 1))
print(solve(-x -2 * y >= -1))
print(solve(1 >= x + 2 * y))
print(solve(1 - x >= 2 * y))

上は$x+2y \le 1$を様々な形で表した場合にどうなるかを調べていて、結果は下のようになる:

2.0  :  x + 2y <= 1
-2.0 : -x - 2y >= -1
2.0  :       1 >= x + 2y
-2.0 :   1 - x >= 2y

シャドウプライスを右辺の増分に対して定義されると考える場合、上から2.0, -2.0, -2.0, -2.0なので、3つ目だけ符号が逆ということになる。これだけが異なるというのは自然な説明が難しく、驚きを小さくするという原則に反しているようにも思える。PuLPのLpAffineExpressionの実装を見るとこうなっている理由がわかる:

class LpAffineExpression(_DICT_TYPE):
    ...
    def __le__(self, other):
        return LpConstraint(self - other, const.LpConstraintLE)

    def __ge__(self, other):
        return LpConstraint(self - other, const.LpConstraintGE)

__le__は<=の実装、__ge__は>=の実装である。a <= bはa - b <= 0と等価であると読める。これはシャドウプライスの観点でも等価である。しかし、aがfloatの場合、Pythonは最初にfloatの__le__を調べ、floatにはLpAffineExpressionとの比較が実装されていないので次に対称律に従ってbの__ge__を試すというステップを踏む。その結果b - a >= 0と解釈されてしまう。上の例でも、1 >= x + 2yの代わりにpl.LpAffineExpression(1) >= x + 2yとあらかじめ変換しておけばシャドウプライスが意図通りの符号になる。

つまりこの問題は、Pythonが比較演算子を対称律の成り立つ述語として想定しているところを、述語でもない上にa >= bとb <= aで意味が異なりうるような使い方をしたことが原因で生じていると言え、解決が難しい。実際、OR-ToolsやPyomoなど他のモデリングツールも同じ問題を抱えているようだ。

OR-Tools:

from collections.abc import Callable
from ortools.linear_solver.pywraplp import Solver, Constraint, Variable


def solve(constr_fn: Callable[[Variable, Variable], Constraint]) -> float:
    solver = Solver.CreateSolver("GLOP")
    x = solver.NumVar(0, Solver.Infinity(), "x")
    y = solver.NumVar(0, Solver.Infinity(), "y")
    solver.Maximize(2 * x + y)
    constr = constr_fn(x, y)
    solver.Add(constr, "constr")
    solver.Solve()
    solver.ComputeConstraintActivities()
    return next(c.dual_value() for c in solver.constraints() if c.name() == "constr")


print(solve(lambda x, y: x + 2 * y <= 1))
print(solve(lambda x, y: -x - 2 * y >= -1))
print(solve(lambda x, y: 1 >= x + 2 * y))
print(solve(lambda x, y: 1 - x >= 2 * y))

# 2.0
# -2.0
# 2.0
# -2.0

Pyomo:

# Pyomo
from collections.abc import Callable
import pyomo.environ as pyo


def solve(constr_fn: Callable[[pyo.ScalarVar, pyo.ScalarVar], pyo.Constraint]) -> float:
    model = pyo.ConcreteModel()
    model.x = pyo.ScalarVar(bounds=(0, None))
    model.y = pyo.Var(bounds=(0, None))
    model.obj = pyo.Objective(expr=2 * model.x + model.y, sense=pyo.maximize)
    model.constr = pyo.Constraint(expr=constr_fn(model.x, model.y), name="constr")
    model.dual = pyo.Suffix(direction=pyo.Suffix.IMPORT)
    pyo.SolverFactory("glpk").solve(model)
    return model.dual[model.constr]


print(solve(lambda x, y: x + 2 * y <= 1))
print(solve(lambda x, y: -x - 2 * y >= -1))
print(solve(lambda x, y: 1 >= x + 2 * y))
print(solve(lambda x, y: 1 - x >= 2 * y))

# 2.0
# -2.0
# 2.0
# 2.0

OR-ToolsはPuLPと同じだが、Pyomoはさらに$1-x \ge 2y$の場合も予想と違っている。

結局、意図通りの符号のシャドウプライスが欲しい場合は使っているツールの実装の詳細に気を配る必要があるということになる。上で見たようにPythonの基本的な挙動に関わる問題なので、完全に「自然」に動くツールがあるとは思えない。

ちなみにAMPLでも試すと、2, -2, -2, -2となって意図通りになる:

option solver 'gurobi';

var x >= 0;
var y >= 0;

maximize obj: 2*x + y;

subject to c1: x + 2*y <= 1;

solve;
display c1.dual;
drop c1;

subject to c2: -x - 2*y >= -1;
solve;
display c2.dual;
drop c2;

subject to c3: 1 >= x + 2*y;
solve;
display c3.dual;
drop c3;

subject to c4: 1 - x >= 2*y;
solve;
display c4.dual;
drop c4;

  1. シャドウプライスはWikipediaでは経済学の言葉として説明されている。感度 (sensitivity) という用語もあるが、これは言葉の指し示す範囲が広く、制約や目的関数の係数を増減させた場合の効果もそう呼ばれていると認識している。最適化系の文献ではシャドウプライスは感度の一種として特別な名前がついていないこともある。また、後述するように、ラグランジュ双対問題との強双対性がある場合はシャドウプライスは双対変数の値に一致するので、dual multiplierとかLagrange multiplierといった概念と同一視されていることも多いと思う。 ↩︎

  2. 例えばMartin, R. Kipp. “Large scale linear and integer optimization.” のCorollary 2.49 ↩︎

  3. Boyd, Stephen P. and Lieven Vandenberghe. “Convex Optimization.” の 5.6. Perturbation and sensitivity analysis などを参照 ↩︎

  4. ただし、ドキュメントには説明がない。関連するPR (#202, #781) やイシューを眺める限り、半公式の機能とみなされているとは思う。 ↩︎

投稿者