双対変数を直感的に理解したい

この記事は 数理最適化Advent Calendar 2023 の9日目の記事です。

導入

数理最適化、特に線形計画問題において、どの教科書にも載っている重要なトピックの一つに双対問題、および双対変数がある。

元の最適化問題の最適値の良い上界、および下界を求める手法の一つとして導入され、数々のアルゴリズムにおいて重要な概念ではあるものの、しばしば形式的に導出されるため双対変数そのものが持つ意味について直感的に分かりづらい点もある。

本記事では、数学的に厳密な話をせずに、双対変数、および種々の性質を直感的に捉えてみることを考える。

双対問題と双対変数

まず、双対問題と双対変数について簡単にまとめる。各定理の詳細や導出に関しては、文献^[1]を始めとした様々な教科書が参考になる。

次のような不等式制約を含む最小化問題を考える。

$$\begin{align} & \min_{x} \sum_{j=1}^{n} c_j x_j \\ \mathrm{s.t. }\ & \sum_{j=1}^{n} a_{kj}x_j \leq b_k \ \forall k\ \\ & x_j \geq 0 \ \forall j \end{align}$$

この時、双対問題は次のような最大化問題の形となる。

$$\begin{align} & \max_{y} \sum_{k=1}^{m} b_k y_k \\ \mathrm{s.t. }\ & \sum_{k=1}^{m} a_{kj}y_k \geq c_j \ \forall j \\ & y_k \leq 0 \ \forall k \end{align}$$

ここで、変数 $y_k$が双対変数である。

元の問題 (主問題) と双対問題の間には以下の性質があることが知られている。

• 双対変数は主問題の各制約 (インデックス $k$)に対応しており、逆に主問題の変数は双対問題の各制約 (インデックス $j$)に対応する。

• 主問題の制約が不等式である場合、対応する双対変数に$y_k \leq 0$ のような非負、または非正の制約がつく。

• 主問題、双対問題のいずれかが実行不能 (実行可能解が存在しない)であるとき、他方の問題は非有界 (解が$\pm \infty$)となる。

• 主問題、双対問題に実行可能解が存在する時、それらの解 $x^{*}, y^{*}$ に対して以下の不等式が成り立つ (弱双対定理)

  $$\sum_{j=1}^{n}c_jx^{*}_{j} \geq \sum_{k=1}^{m}b_ky^{*}_{k}$$

  • 特に線形計画問題では、それらの解 $x^{*}, y^{*}$ が最適値である場合は、この不等式は等式となる (強双対定理)

    $$\sum_{j=1}^{n}c_jx^{*}_{j} = \sum_{k=1}^{m}b_ky^{*}_{k}$$

このため、双対問題は以下のような用途に使われる。

• 弱双対定理により元の最小化問題の下限を見積もることができるため、得られた解がどこまで最適解に近いかを調べることができる。
• 主問題の変数が極めて多く、制約条件が比較的少ない場合に、双対問題に移行することで制約が多く、変数が比較的少ない問題に書き換えることができるため、このテクニックを用いて求解を容易にすることができるケースがある。(e.g. 列生成法)

少し実験してみる

ここで、数理最適化のモデラであるPythonMIPを使って試しに次の問題を解き、双対変数を取得してみる。

ここで、各制約のコロンの後の変数 $\lambda_1, \lambda_2$ はそれぞれの制約に対応する双対変数を表す。

$$\begin{align} & \min_{x,y} & 3x &+ 4y \\ &\mathrm{s.t. }\ & 2x &+ y & \geq & 20 &: \lambda_1 \\ & & x&+2y & \geq & 30 &: \lambda_2 \\ & & x&,y &\geq &0 \ \end{align}$$

PythonMIPは以下のようにpip経由でインストールでき、上の問題を解くコードは以下である。

pip install mip -U

from mip import Model, xsum, MINIMIZE, CONTINUOUS

# モデルの作成
model = Model(sense=MINIMIZE)

# 変数の定義（非負）
# 連続変数として定義しないと双対変数は取得できない。
x = model.add_var(var_type=CONTINUOUS, name="x")
y = model.add_var(var_type=CONTINUOUS, name="y")

# 目的関数の定義：最大化 3x + 4y
model.objective = 3 * x + 4 * y

# 制約の追加
model += (2 * x + y >= 20, "const1")
model += (x + 2 * y >= 30, "const2")

# 問題を解く
model.optimize()

# 結果の取得（最適解と目的関数の値）
optimal_x = x.x
optimal_y = y.x
optimal_value = model.objective_value

print("x = ", optimal_x, ", y = ", optimal_y, ", value = ", optimal_value)
# get dual variables
print("const1 dual variable: ", model.constr_by_name("const1").pi)
print("const2 dual variable: ", model.constr_by_name("const2").pi)

結果は次のようになる。

x =  3.333333333333333 , y =  13.333333333333334 , value =  63.333333333333336
const1 dual variable:  0.6666666666666667
const2 dual variable:  1.6666666666666665

この問題では主問題の最適解が$x=3.33, y = 13.33$ で最適値が63.33となり、双対問題の解が$\lambda_1=0.67, \lambda_2=1.67$となることがわかる。

ここで、元の問題の制約の右辺を微小に (+0.1)変化させて再び問題を解いてみよう。

$$\begin{align} & \min_{x,y} & 3x &+ 4y \\ &\mathrm{s.t. }\ & 2x &+ y & \geq & 20+0.1 &: \lambda_1 \\ & & x&+2y & \geq & 30 &: \lambda_2 \\ & & x&,y &\geq &0 \ \end{align}$$

上のコードを書き換え、再度実行すると最適値が63.4と、先程の値に比べ0.067だけ変化していることがわかる。

x =  3.400000000000001 , y =  13.299999999999999 , value =  63.4
const1 dual variable:  0.6666666666666667
const2 dual variable:  1.6666666666666665

また、変化量の比率が、

$$\frac{最適値の変化量}{制約式の右辺の変化量} = \frac{0.067}{0.1} = 0.67 = \lambda_1$$

と、双対変数と全く同じ値になっていることもわかるだろう。

もう一方の制約条件も同様の計算をすると$\lambda_2$ と同じになることがすぐわかる。

双対変数と変化量の関係

上の実験結果を一般化して考えてみる。

まず、元の問題の制約式の右辺を$\delta b_k$ だけずらした次の問題を考えてみよう。

$$\begin{align} & \min_{x} \sum_{j=1}^{n} c_j x_j \\ \mathrm{s.t. }\ & \sum_{j=1}^{n} a_{kj}x_j \leq b_k + \delta b_k \ \forall k\ \\ & x_j \geq 0 \ \forall j \end{align}$$

この問題の双対問題を考えると、$b_k$を$b_k + \delta b_k$に置き換えるだけであるため、以下のようになる。

$$\begin{align} & \max_{y} \sum_{k=1}^{m} (b_k + \delta b_k) y_k \\ \mathrm{s.t. }\ & \sum_{k=1}^{m} a_{kj}y_k \geq c_j \ \forall j \\ & y_k \leq 0 \ \forall k \end{align}$$

これにより、元のコスト関数の最適値を$z$とすると、強双対定理より

$$\sum_{j=1}^{n}c_jx^{*}_{j} = \sum_{k=1}^{m}b_ky^{*}_{k}$$

であるため、変化分$\delta z$ は以下のようになる。

$$\delta z = \sum_{k=1}^{m}y_k \delta b_k$$

すなわち、双対変数はある制約式の右辺を微小に変化させた際の、対応するコスト関数の変化分の比率を表す量となることがわかる。

特に、ある一つだけの制約の右辺を$\delta b_k$ だけ微小に変化させた際、コスト関数の最適値が$\delta z$ だけ変化したとすると、対応する双対変数は

$$y_k = \frac{\delta z}{\delta b_k}$$

により得ることができる。

この性質自体は感度分析において用いられ、対象となる制約条件が元のコスト関数にどれだけ重要な影響を及ぼすかを測定する指標としてよく利用されている。^[2]

以下の節では、双対変数に関する種々の性質を上の概念を用いて見ていこうと思う。

不等式制約に対応する双対変数

主問題に不等式制約が含まれる際、双対問題において対応する双対変数に非正制約、あるいは非負制約がかかる。この性質を見ていく。

今、制約に$\leq$の不等号が含まれる次の問題を考える。

$$\begin{align} & \min_{x} \sum_{j=1}^{n} c_j x_j \\ \mathrm{s.t. }\ & \sum_{j=1}^{n} a_{kj}x_j \leq b_k \ \forall k\ \\ & x_j \geq 0 \ \forall j \end{align}$$

この制約式にて、制約式の右辺を$\delta b_k (> 0)$ だけ変化させると

$$\sum_{j=1}^{n} a_{kj}x_j \leq b_k + \delta b_k$$

この変化は元あった制約を緩める方向へと変化させていることがわかる。そのため主問題の決定変数の取りうる自由度は大きくなり、コスト関数はより小さい値 ($\delta z \leq 0$)を取れるようになる。

$\delta b_k >0, \delta z \leq 0$ であるため、対応する双対変数の符号は

$$y_k = \frac{\delta z}{\delta b_k} \leq 0$$

となるため、双対変数は必ず0以下の値となることがわかる。

同様に、制約に$\geq$の不等号が含まれる場合、

$$\sum_{j=1}^{n} a_{kj}x_j \geq b_k + \delta b_k$$

の変化は元あった制約をきつくする方向へと変化させていることがわかる。そのため主問題の決定変数の取りうる自由度は小さくなり、コスト関数はより大きい値 ($\delta z \geq 0$)を取らざるをえない。

この場合、$\delta b_k >0, \delta z \geq 0$ であるため、

$$y_k = \frac{\delta z}{\delta b_k} \geq 0$$

となり、双対変数は必ず0以上の値となることがわかる。^[3]

双対変数が0になるケース

ここで、双対変数が0となるケースに考えてみよう。

$$\delta z = \sum_{k=1}^{m}y_k \delta b_k$$

において、双対変数$y_k$が0であれば、対応する主問題の制約の右辺に微小な変化を加えてもコスト関数に影響がないことがわかる。問題から取り除いても最適解に影響しない制約 (=有効でない制約) についてこの性質が生じる。^[4]

双対変数が無限大に発散するケース

主問題が実行不能な問題で双対問題に実行可能解がある時^[5]、その双対問題の解は非有界になる。この性質をみてみよう。

今、次の問題を考える。

$$\begin{align} & \min_{x} \sum_{j=1}^{n} c_j x_j \\ \mathrm{s.t. }\ & \sum_{j=1}^{n} a_{kj}x_j \leq b_k \ \forall k\ \\ & x_j \geq 0 \ \forall j \end{align}$$

この問題を常に実行可能な問題にするため、以下の書き換えを行ってみる。

$$\begin{align} & \min_{x, s} \sum_{j=1}^{n} c_j x_j + M\sum_{k} s_k \\ \mathrm{s.t. }\ & \sum_{j=1}^{n} a_{kj}x_j \leq b_k + s_k \ \forall k\ \\ & x_j \geq 0 \ \forall j \\ & s_k\geq 0 \ \forall k \\ \end{align}$$

ここで、$M$は比較的大きな正の定数で、$s_k$がいわゆるslack変数である。

もし元の問題が実行可能解を持つのであれば、slack変数$s_k$がコスト関数の第2項によりペナルティをかけられているため、$s_k = 0$ となる解が出てくるはずである。実際に

$$\begin{align} & \min_{x,y} & -3x &+ Ms \\ &\mathrm{s.t. }\ &x & \leq 20+s \\&&x,s &\geq 0 \ \end{align}$$

を計算してみるとそのとおりとなる。

from mip import Model, xsum, MINIMIZE, CONTINUOUS

# モデルの作成
model = Model(sense=MINIMIZE)

# 変数の定義（非負）
# 連続変数として定義しないと双対変数は取得できない。
x = model.add_var(var_type=CONTINUOUS, name="x")
# slack変数
s = model.add_var(var_type=CONTINUOUS, name="s")
# 大きい定数
M = 10000
# 目的関数の定義：最小化 -3x + M * s
model.objective = -3 * x + M * s

# 制約の追加
model += (x <= 20 + s, "const1")

# 問題を解く
model.optimize()

# 結果の取得（最適解と目的関数の値）
optimal_x = x.x
optimal_s = s.x
optimal_value = model.objective_value

print("x = ", optimal_x, ", s = ", optimal_s, ", value = ", optimal_value)
# get dual variables
print("const1 dual variable: ", model.constr_by_name("const1").pi)

x =  20.0 , s =  0.0 , value =  -60.0
const1 dual variable:  -3.0

元の問題に実行可能解がない場合、スラック変数は0で無い値を持つ。ここで、

$$a_{kj}x_j \leq b_k + s_k + \delta b_k$$

のように制約式の右辺に微小な変化を加えると、連動してスラック変数が変動し、コスト関数の第2項の$M \sum_{k}s_k$ の値が変化する。$M$は大きな定数なのでコスト関数の変化値、ひいては双対変数が$M$に依存する結果となる。実際に以下の(元の問題が実行不能である)問題で試してみると、双対変数がかなり大きなスケールの値となることがわかる。

$$\begin{align} & \min_{x,y} & -3x &+ Ms \\ &\mathrm{s.t. }\ &x & \leq -5+s \\&&x,s &\geq 0 \ \end{align}$$

from mip import Model, xsum, MINIMIZE, CONTINUOUS

# モデルの作成
model = Model(sense=MINIMIZE)

# 変数の定義（非負）
# 連続変数として定義しないと双対変数は取得できない。
x = model.add_var(var_type=CONTINUOUS, name="x")
# slack変数
s = model.add_var(var_type=CONTINUOUS, name="s")
# 大きい定数
M = 10000
# 目的関数の定義：最小化 -3x + M * s
model.objective = -3 * x + M * s

# 制約の追加
model += (x <= -5 + s, "const1")

# 問題を解く
model.optimize()

# 結果の取得（最適解と目的関数の値）
optimal_x = x.x
optimal_s = s.x
optimal_value = model.objective_value

print("x = ", optimal_x, ", s = ", optimal_s, ", value = ", optimal_value)
# get dual variables
print("const1 dual variable: ", model.constr_by_name("const1").pi)

x =  0.0 , s =  5.0 , value =  50000.0
const1 dual variable:  -10000.0

元の問題はこのスラック変数のある問題にて、$s_k = 0$ という条件を加えたものに等しいので、$M$を極めて大きい値にし、$s_k$を必ず0以外の値にしないようなペナルティをかけている状態と見ることができる。この時、$M$を大きくすればするほど対応する双対変数のスケールも際限なく大きくなることが見て取れる。

これにより元の問題が実行不能である場合、双対変数が非有界になる現象を直感的に見ることができたと思う。

双対変数を使ったアルゴリズム

最後に、この双対変数を使って少し発展的なトピックを扱ってみる。

今、下のような構造を持った問題を扱ってみよう。

$$\begin{align} & \min_{x,y} \sum_{j=1}^{n} c_j x_j \\ \mathrm{s.t. }\ & \sum_{j=1}^{n} a_{kj}x_j + y \leq b_k \ \forall k\ \\ & x_j \geq 0 \ \forall j \end{align}$$

決定変数に$y$が一部含まれているが、もしこの$y$がなかったら解くのが簡単な問題なのになぁ、といったシチュエーションは数理最適化にて頻繁に起こる。例えば

• $y$がなければ残った問題が複数の小さな問題に分解できる
• $y$が離散変数 (バイナリ変数等)であるため、$y$さえ固定できてしまえば残りの問題が線形計画問題として解ける

といった状況である。そのため、$x$と$y$を別々に解きたい。

一つ思いつく方法としては、とりあえず$y$をいろいろ動かしてみて

$$\begin{align} & g(y) = \min_{x} \sum_{j=1}^{n} c_j x_j \\ \mathrm{s.t. }\ & \sum_{j=1}^{n} a_{kj}x_j + y \leq b_k \ \forall k\ \\ & x_j \geq 0 \ \forall j \end{align}$$

で定義される$y$に対する関数$g(y)$ の最小値 $\min_{y} g(y)$ を見つけることである。しかし、(特に$y$が多次元変数等の場合) 各々の$y$に対して関数 $g(y)$ の値を一つ求めるのに最適化問題を解かないといけないため、正確に復元することはかなり時間がかかってしまう。

ここで双対変数が「ある制約式の右辺を微小に変化させた際の、対応するコスト関数の変化分の比率」であることを思い出すと、もし

$$y = y^{*}$$

のように$y$を固定し、これを上の問題の制約に入れた場合、この制約に対する双対変数は$y = y^{*} + \delta y$ と動かした際の関数 $g(y)$ の変化値となるはずだから、$g(y)$ の$y = y^{*}$ における傾きを指す値となることがわかる。

すなわち次の問題を解くことにより、得られた$y = y^{*}$に対する双対変数がその点における$g(y)$ の傾きとなる。

$$\begin{align} & \min_{x,y} \sum_{j=1}^{n} c_j x_j \\ \mathrm{s.t. }\ & \sum_{j=1}^{n} a_{kj}x_j + y \leq b_k \ \forall k\ \\ & y = y^{*}: \lambda_1 \\ & x_j \geq 0 \ \forall j \end{align}$$

このため、様々な場所で$y$ を固定することで、元の関数である$g(y)$ の概形を復元することができる。

さらに、$y$ を闇雲に固定するのではなく、固定する手順を次のように決めることもできる。

1. $y$について次の何の制約もない最適化問題を解く。$y$の値は何か適当な値となる
   変数$\theta$はこの後のステップで用いる。

   $$\begin{align} & \min_{y, \theta} \theta \\ \end{align}$$

2. 1で得られた$y$の値を$y^{*(1)}$とし、次の問題を$y$を固定した状態で解き、双対変数$\lambda_1^{(1)}$、およびコスト関数$v^{(1)}$を得る。

   $$\begin{align} & \min_{x,y} \sum_{j=1}^{n} c_j x_j \\ \mathrm{s.t. }\ & \sum_{j=1}^{n} a_{kj}x_j + y \leq b_k \ \forall k\ \\ & y = y^{*(1)}: \lambda_1^{(1)} \\ & x_j \geq 0 \ \forall j \end{align}$$

3. 1に次の制約を追加し、再び$y$ について解く。

   $$\begin{align} &\min_{y, \theta} \theta \\ \mathrm{s.t. }\ & \theta \geq v^{(1)} + \lambda_1^{(1)}(y - y^{*(1)}) \end{align}$$

4. 3で得られた$y$の値を$y^{*(2)}$とし、2を繰り返す。

これを繰り返すことで、$y$に関する問題に次のように逐次的に制約が追加される形となる。

$$\begin{align} &\min_{y, \theta} \theta \\ \mathrm{s.t. }\ & \theta \geq v^{(1)} + \lambda_1^{(1)}(y - y^{*(1)}) \\ & \theta \geq v^{(2)} + \lambda_1^{(2)}(y - y^{*(2)}) \\ & \theta \geq v^{(3)} + \lambda_1^{(3)}(y - y^{*(3)}) \\ \cdots \\ \end{align}$$

追加された制約により得られる実行可能解空間は、下の図の塗りつぶした領域となるので、この処理を繰り返すとやがて最小値 $\min_{y} g(y)$ を見つけることができる。
これにより元の問題の最適解を求めることができる。

実はこのアルゴリズムはBenders分解法^[6]という名前で知られており、今回のように直接問題を解くのが難しいケースにおいて各変数に対して交互に解く手法として知られている。^[7]

まとめ

この記事では、双対変数が「ある制約式の右辺を微小に変化させた際の、対応するコスト関数の変化分の比率」に相当することを見て、これを用いて双対変数の種々の性質を見てきた。
さらに、双対変数を使うことで数理最適化において重要なアルゴリズムを構築できることも見た。

双対変数について教科書で定義や性質を学んだものの、双対変数の直感的な意味そのものについてあまり注目したことがなかった読者もいたのではないだろうか。

今後双対変数に関するトピックが出てきた際、

• この双対変数は0だから、対応する制約は問題に寄与していなさそうだ。
• 双対変数が大きい値を持つから、この制約がコストにめっちゃ効いてきそう。

のように、実務的なシチュエーションにおいてもより双対変数を身近に感じてもらえれば幸いである。

脚注

1. 梅谷俊治,「しっかり学ぶ数理最適化」, 講談社 (2020)

2. 感度分析の文脈では、制約条件が間接的にコスト関数に与える影響を示すものという意図で双対変数をshadow price と呼ぶこともある。

3. 対象とする最適化問題が最小化ではなく最大化である場合、双対変数の符号の関係はこの記事とは異なる。ただ、この記事と同じ手法で確かめることができる。

4. 詳しくは相補性定理を参照

5. 主問題、双対問題とともに実行不能であるケースも存在することに注意。

6. Benders, J. F. (Sept. 1962), "Partitioning procedures for solving mixed-variables programming problems", Numerische Mathematik 4(3): 238–252.

7. このformulationは"Decomposition Techniques in Mathematical Programming (Springer-Verlag, 2006)."を参考にしている。