2013-04-09 [長年日記]
λ. 強双対性を使って線形計画を制約充足問題に還元する
Linear and Discrete Optmization の双対性に関する部分(参考: 第五週目メモ)を復習していて気になったのだけど、強双対性の条件 cT x* = bT y* って線形式で書けてしまっているのよね。 ということは主問題と双対問題の制約条件と強双対性の連言をとっても線形(QF_LRA)の範囲に収まっていて、これの解(実行可能解)を求めれば、目的関数に関する最適化をしなくても最適解が求まることになるはず。 ということで試してみた。
主問題
簡単な問題として、Linear and Discrete Optmization で最初の例として使われていた、以下の例を考える。
MAXIMIZE 100 x1 + 125 x2 SUBJECT TO 3 x1 + 6 x2 <= 30 8 x1 + 4 x2 <= 44 BOUNDS 0 <= x1 <= 5 0 <= x2 <= 4 END
最適解は (x1,x2)=(4,3) で、目的関数値は775。
双対問題
この問題の双対問題は以下のようになる。
MINIMIZE 30 y1 + 44 y2 + 5 y3 + 4 y4 SUBJECT TO 3 y1 + 8 y2 + y3 >= 100 6 y1 + 4 y2 + y4 >= 125 BOUNDS 0 <= y1 0 <= y2 0 <= y3 0 <= y4 END
最適解は (y1,y2,y3,y4) = (50/3, 25/4, 0, 0) で、目的関数値は775。
両方を合わせたもの
MAXIMIZE 0 x1 SUBJECT TO 3 x1 + 6 x2 <= 30 8 x1 + 4 x2 <= 44 3 y1 + 8 y2 + y3 >= 100 6 y1 + 4 y2 + y4 >= 125 100 x1 + 125 x2 - 30 y1 - 44 y2 - 5 y3 - 4 y4 = 0 BOUNDS 0 <= x1 <= 5 0 <= x2 <= 4 0 <= y1 0 <= y2 0 <= y3 0 <= y4 END
(x1,x2,y1,y2,y3,y4) = (4, 3, 50/3, 25/4, 0, 0) はこの問題の(実行可能)解で、このときの目的関数値はどちらも775。
Z3の場合
最適化機能を持たない充足可能性ソルバの例として、SMTソルバのZ3で試してみる（オンラインでも Z3 online at RiSE4Fun で試せる）。
(declare-fun x1 () Real) (declare-fun x2 () Real) (define-fun obj-primal () Real (+ (* 100.0 x1) (* 125.0 x2))) (assert (<= (+ (* 3.0 x1) (* 6.0 x2)) 30.0)) (assert (<= (+ (* 8.0 x1) (* 4.0 x2)) 44.0)) (assert (and (<= 0.0 x1) (<= x1 5.0))) (assert (and (<= 0.0 x2) (<= x2 4.0))) (declare-fun y1 () Real) (declare-fun y2 () Real) (declare-fun y3 () Real) (declare-fun y4 () Real) (define-fun obj-dual () Real (+ (* 30.0 y1) (* 44.0 y2) (* 5.0 y3) (* 4.0 y4))) (assert (>= (+ (* 3.0 y1) (* 8.0 y2) y3) 100.0)) (assert (>= (+ (* 6.0 y1) (* 4.0 y2) y4) 125.0)) (assert (<= 0.0 y1)) (assert (<= 0.0 y2)) (assert (<= 0.0 y3)) (assert (<= 0.0 y4)) (assert (= obj-primal obj-dual)) (check-sat) ; => sat (get-value (obj-primal x1 x2)) ; => ((obj-primal 775.0) (x1 4.0) (x2 3.0)) (get-value (obj-dual y1 y2 y3 y4)) ; => ((obj-dual 775.0) (y1 (/ 50.0 3.0)) ; (y2 (/ 25.0 4.0)) (y3 0.0) (y4 0.0))
感想とか
整数上の場合だと「解が一つ見つかった後に、目的関数値がそれよりも真に良いということを表す制約を追加して、再度解き直す」ということを繰り返すことで、有界な場合には最適化を行えるのに対して、実数上だと有界であってもZeno状態みたいなことに陥って、最適解に到達できない場合があるはず。 なので、実数の場合には充足可能性ソルバでは最適化問題は解けないと思っていたけれど、こういうやり方をすれば解けるんだなぁ。 ということで面白かった。
まあ、浮動小数点数で計算するような場合には、誤差蓄積で解が見つからなくなるような場合もあるような気はするけれど。
2013-04-04 [長年日記]
λ. ゴモリーの小数カットの導出メモ
Outline of an algorithm for integer solutions to linear programs でのゴモリーカット(Gomory's fractional cut)が、割と天下りに定義して、そのあと必要な性質を満たすことを証明しているけれど、これがちょっと気に入らなかったので、演繹的に導出してみた。(An algorithm for the mixed integer problem のゴモリーの混合整数カットの方ではないので注意)
準備
N を(0を含む)自然数の集合とし、(R, ≤R), (Z, ≤Z) をそれぞれ実数と整数の順序集合とする。 また、 I : (Z, ≤Z) → (R, ≤R) を順序集合間の通常の埋め込みとするが、たいてい省略する。
⌈-⌉, ⌊-⌋ : (R, ≤R) → (Z, ≤Z) をそれぞれ正の無限大方向への切り上げ(ceiling)と負の無限大方向への切り捨て(floor)とする。
随伴、モナド、コモナド、双対性
⌈-⌉ と ⌊-⌋ はそれぞれ I の左随伴と右随伴になっている。すなわち、以下が成り立つ。
- ∀a∈R, b∈Z. ⌈a⌉ ≤Z b ⇔ a ≤R I(b)
- ∀a∈Z, b∈R. I(a) ≤R b ⇔ a ≤Z ⌊b⌋
したがって、 I∘⌈-⌉, I∘⌊-⌋ : (R, ≤R) → (R, ≤R) はそれぞれ (R, ≤R) 上のモナド、コモナドとなっており、以下が成り立つ。
- ⌊a⌋ ≤ a ≤ ⌈a⌉
- ⌊⌊a⌋⌋ = ⌊a⌋
- ⌈a⌉ = ⌈⌈a⌉⌉
また、 ⌈-⌉∘I, ⌊-⌋∘I : (Z, ≤Z) → (Z, ≤Z) はそれぞれ (Z, ≤Z) 上のコモナド、モナドとなっており、こっちは更に強く、以下が成り立っている。
- ⌈I(a)⌉ = a for all a∈Z
- a = ⌊I(a)⌋ for all a∈Z
また、(-1)×- : (R, ≤R) → (R,≥R) によって、⌈-⌉と⌊-⌋は双対になっており、 ⌈a⌉ = - ⌊- a⌋ が成り立つ。
⌈-⌉ と ⌊-⌋ の分配性
a≤⌈a⌉ および b≤⌈b⌉ より a+b ≤R ⌈a⌉+⌈b⌉ = I(⌈a⌉+⌈b⌉) で、随伴性より ⌈a+b⌉ ≤Z ⌈a⌉+⌈b⌉ 。 双対性より ⌊a+b⌋ ≥Z ⌊a⌋+⌊b⌋ 。
さらに、帰納法により ∀n∈N. ⌈n×a⌉ ≤Z n×⌈a⌉ もいえる。 同じく、双対性より ∀n∈N. ⌊n×a⌋ ≥Z n× ⌊a⌋
ゴモリーの小数カットの導出
xi, y ∈ N が
- Σi ai xi + y = b …… (1)
を満たすとする。すると、
- Σi ai xi + y ≥ b … (2)
- Σi ai xi + y ≤ b … (3)
(3)の両辺に⌊-⌋を適用すると、⌊-⌋ の単調性より
- ⌊Σi ai xi + y⌋ ≤ ⌊b⌋ … (4)
⌊-⌋ の分配性と ⌊y⌋=y より
- ⌊Σi ai xi + y⌋ ≥ Σi ⌊ai xi⌋ + ⌊y⌋ ≥ Σi ⌊ai⌋ xi + ⌊y⌋ = Σi ⌊ai⌋ xi + y … (5)
(4), (5)に推移律を適用して
- Σi ⌊ai⌋ xi + y ≤ ⌊b⌋ … (6)
両辺に -1 をかけて
- Σi (- ⌊ai⌋) xi - y ≥ -⌊b⌋ … (7)
(2) と (7) の両辺を足し合わせて、
- Σi (ai - ⌊ai⌋) xi ≥ b - ⌊b⌋
■
2013-04-02 [長年日記]
λ. 「異様に難しいドラクエの謎解き」を整数線形計画問題として解く
「異様に難しいドラクエの謎解き」をモデル検査器で解く の続き。 元の「人生を振り返る - 異様に難しいドラクエの謎解き」に「押す順番はどうでも良くて、四隅のボタンを押す回数の問題でしたね」とあって、この洞察を使えば整数線形計画問題として定式化して、もっとスマートに解けるなぁ、と思った。
整数線形計画問題による定式化
x1, x2, x3, x4 をそれぞれ左上、右上、右下、左下の石を押す回数とすると、最短の手順を求める問題は、以下の最適化問題として定式化できる。
Minimize
x1 + x2 + x3 + x4
Subject To
3 + x1 + x2 + x3 ≡ 0 mod 4 (左の像の条件)
3 + x2 + x3 + x4 ≡ 0 mod 4 (上の像の条件)
x3 + x4 + x1 ≡ 0 mod 4 (右の像の条件)
x4 + x1 + x2 ≡ 0 mod 4 (下の像の条件)
x1, x2, x3, x4 ≧ 0
x1, x2, x3, x4 ∈ Z
追加的な変数を導入してmodを取り除けば、これは整数線形計画問題になる。その結果の問題をLPファイル形式にすると以下ようになる(baroque_tower.lp)。
Minimize x1 + x2 + x3 + x4 Subject To c1: x1 + x2 + x3 - 4 n1 = -3 c2: x2 + x3 + x4 - 4 n2 = -3 c3: x3 + x4 + x1 - 4 n3 = 0 c4: x4 + x1 + x2 - 4 n4 = 0 Bounds 0 <= x1 0 <= x2 0 <= x3 0 <= x4 0 <= n1 0 <= n2 0 <= n3 0 <= n4 General x1 x2 x3 x4 n1 n2 n3 n4 End
MIPソルバで解いてみる
これをMIPソルバのlpsolveを使って解かせてみる。
> lp_solve.exe -rxli xli_CPLEX.dll baroque_tower.lp set_XLI: Successfully loaded 'xli_CPLEX.dll' Value of objective function: 6.00000000 Actual values of the variables: x1 1 x2 2 x3 2 x4 1 n1 2 n2 2 n3 1 n4 1
ということで、(x1,x2,x3,x4)=(1,2,2,1)という最適解が得られ、これは左上を1回、右上と右下を2回ずつ、左下を一回という、前回得られたものと同じ解が得られる。
実は最初はGLPKでも解こうと思ったけれど、GLPKだと、有界でない整数変数のせいで無限にブランチして止まらないっぽい。各変数に適当な上界を加えれば解ける。 あと、自分が以前に書いたオモチャのMIPソルバでも楽勝で解けた。
最適解の一意性
最適解の一意性については、MIPソルバを使うよりも、論理結合子を自由に使えるSMTソルバの方が便利なので、SMTソルバのZ3で確認してみることにする。
以下の内容(baroque_tower.smt2)をZ3に与える（オンラインでも Z3 online at RiSE4Fun で試せる）ことで、6手の他の解が存在しないことを確認できる。 詳しくは以下のコメント中にて。
; baroque_tower.smt2 (set-option :produce-models true) (set-logic QF_LIA) ; 変数の宣言 (declare-fun x1 () Int) (declare-fun x2 () Int) (declare-fun x3 () Int) (declare-fun x4 () Int) ; 前述の制約条件の追加 (assert (>= x1 0)) (assert (>= x2 0)) (assert (>= x3 0)) (assert (>= x4 0)) (assert (= (mod (+ 3 x1 x2 x3) 4) 0)) (assert (= (mod (+ 3 x2 x3 x4) 4) 0)) (assert (= (mod (+ x3 x4 x1) 4) 0)) (assert (= (mod (+ x4 x1 x2) 4) 0)) ;; まずは適当に解を求めてみる (check-sat) ; sat と出力され、解が存在することが分かるので、解を表示させる (get-value (x1 x2 x3 x4)) ; ((x1 1) ; (x2 4946) ; (x3 4990) ; (x4 9749)) ;; 既に最適解での目的関数値が6と分かっているので、その条件を加えて再度解いてみる (assert (= (+ x1 x2 x3 x4) 6)) (check-sat) ; 今度もsatと出力され、解が存在することが分かるので、 解を表示させる (get-value (x1 x2 x3 x4)) ; ((x1 1) ; (x2 2) ; (x3 2) ; (x4 1)) ;; 更に、今見つかった最適解を除外する制約を加えて、再度解いてみる (assert (not (and (= x1 1) (= x2 2) (= x3 2) (= x4 1)))) (check-sat) ; 今度はunsatと出力され、従って今除外された以外の解がないことが分かる
関連リンク
ひょっとしたら関係するかもしれない書籍など
2013-03-27 [長年日記]
λ. 「異様に難しいドラクエの謎解き」をモデル検査器で解く
ドラクエ７のバロックの塔 (バロックタワー) の謎解き (ギミック) が異様に難しい。一応自力で解いてみたけど、手数が結構掛かってしまったので、最適なパターンが無いのか気になった。
この手の問題は、昔書いた「油売り算」と同じで、モデル検査器(Model Checker)を使うと簡単に解ける。
モデル検査器というのは何かというと、並行システムとかのような複雑な状態遷移を持つようなシステムに対して、エラー状態に到達しないだとか、デッドロックしないだとかの性質を検証するツールで、逆にそういうことが起こりうる場合には、それに至る状態遷移列を「反例」としして出力するような機能を持っている。 そこで、この手の問題の場合、状態と可能な手からなる状態遷移系を考えて、ゴールとなるような状態を「エラー」として、モデル検査器にかけてあげると、ゴールに至る状態遷移列が「反例」として得られる。
ということで、モデル検査器SPINを使って6手が最短であることを確認してみる。 まず、以下が検査用の状態遷移系の記述(BaroqueTower.pml)。 非決定的に状態遷移を繰り返し、目的の状態になったら表明違反を発生させるように書いておくのがポイント。
/* 0:左の像, 1:上の像, 2:右の像, 3:下の像 */
byte dir[4]
/* 像を時計回りに90度回転させる */
inline rotate(n)
{
dir[n] = (dir[n] + 1) % 4
}
/* 石から反時計回りに3つの像を回転させる */
/* 0:左上の石, 1:右上の石, 2:右下の石, 3:左下石 */
inline push_button(idx)
{
atomic{
rotate(idx)
rotate((idx+3) % 4)
rotate((idx+2) % 4)
}
}
active proctype main()
{
atomic {
dir[0] = 3; /* 左の像 */
dir[1] = 3; /* 上の像 */
dir[2] = 0; /* 右の像 */
dir[3] = 0; /* 下の像 */
printf("[%d,%d,%d,%d]\n", dir[0], dir[1], dir[2], dir[3]);
}
/* 非決定的に色々試す */
do
:: atomic {
int k;
if
:: true -> k = 0
:: true -> k = 1
:: true -> k = 2
:: true -> k = 3
fi;
push_button(k);
printf("%d -> [%d,%d,%d,%d]\n", k, dir[0], dir[1], dir[2], dir[3]);
}
/* 解けたら表明違反が発生するように */
assert(! (!dir[0] && !dir[1] && !dir[2] && !dir[3]) );
od
}
これを元に検査器を生成してコンパイルする。
% spin -a BaroqueTower.pml % gcc -O2 -DREACH -o pan pan.c
こうして生成された検査器には「-i search for shortest path to error」というオプションがあるので、「pan -i」としてエラーに至る最短パスを探させる。結果は BaroqueTower.pml.trail というファイルに書き込まれる。
% pan -i
得られたパスの情報を元にシミュレーションさせる。
% spin -k BaroqueTower.pml.trail BaroqueTower.pml
[3,3,0,0]
1 -> [0,0,0,1]
1 -> [1,1,0,2]
2 -> [2,2,1,2]
2 -> [3,3,2,2]
3 -> [3,0,3,3]
0 -> [0,0,0,0]
spin: BaroqueTower.pml:46, Error: assertion violated
spin: text of failed assertion: assert(!((((!(dir[0])&&!(dir[1]))&&!(dir[2]))&&!(dir[3]))))
spin: trail ends after 37 steps
#processes: 1
dir[0] = 0
dir[1] = 0
dir[2] = 0
dir[3] = 0
37: proc 0 (main) BaroqueTower.pml:32 (state 29)
1 process created
ということで、 最短パス(の一つ)は「右上→右上→右下→右下→左下→左上」(1,1,2,2,3,0)で、6手未満のパスはないことが確認できる。
まあ、これくらいの規模の問題なら、モデル検査器を使わずに、自分で適当にプログラムを書いて状態遷移系を探索しても、大したことはないけれど。
参考図書など
2013-03-15 [長年日記]
λ. TAPLの日本語訳出版します
すっかり書きそびれていたのですが、Benjamin C. Pierce の Types and Programming Languages (TAPL) の翻訳(共訳)をここしばらくしていたのですが、それがようやくこの3月に出版されることになりました。
しかもトークイベントなるものがあるので、何を喋るのか謎だけれど(^^;、興味のある方がいましたら、是非ご参加ください。
記憶が正しければ、原著について知ったのは、出版されたときに Benjamin C. Pierce のメールを向井先生に見せてもらって知ったのだったと思う。 当時はまさかそれを自分が翻訳することになるとは思っておらず、今更ながら驚いている*1。
ただ、当時はTAPLについて特に興味を持つこともなかったし、その後、型理論系の論文を読むようになった後も、何かの機会に買ってはみたものの、結局積読になってしまっていた。TAPLに書いてあるようなことはある程度は知っているつもりだったので、あまり読む必要性を感じていなかったというのもある。が、今回自分たちで翻訳してみて、これはすごく勿体ないことをしてたなぁ、と思った。 今となっては当然知っていることが多いのだけれど、知っていることについてもすごくしっかり書いてあるし、知らないこともそれなりにあった。 TAPL一冊読めば相当力が付くだろうし、TAPLを読まずにいたことで、昔の自分はずいぶん遠回りをしていたなぁと。
*1 原著は2002年2月の出版だから、当時の自分は学部1年だったわけで、それから10年以上経っていると思うと、ずいぶん遠くまで来たなぁと思う。
