Day 131-Q5 | Seidelのランダム化増分法（2次元線形計画法）

問題

平面上に$N$個の半平面制約$a_ix+b_iy\le c_i$が与えられる。これらすべてを満たす領域の中で、目的関数$p\,x+q\,y$を最小化する点を求めたい。実行可能領域が空ならINFEASIBLE、非有界ならUNBOUNDED、それ以外は最小値（誤差$10^{-6}$以下許容）を出力せよ。

入力形式

N
p q
a_1 b_1 c_1
a_2 b_2 c_2
...
a_N b_N c_N

制約

$1 \le N \le 200$

$-10^6 \le a_i,b_i,c_i,p,q \le 10^6$

$(a_i,b_i) \neq (0,0)$

入出力例

入力例1

出力例1

2.000000

x≥1,y≥1,x+y≤10,x≤5,y≤5 のもとでx+yを最小化 → (1,1)で値2。

入力例2

出力例2

INFEASIBLE

入力例3

1
1 0
0 -1 -1

出力例3

UNBOUNDED

概念図

ヒント（段階的開示）

ヒント1（方向性）

2次元LPは単体法でも解けるが、「LP-type問題」という一般的な枠組みで捉えると乱択増分法により期待O(N)というアルゴリズムが存在する。これがSeidelのアルゴリズムで、Welzlの最小包含円と同じLP-type理論に属する親戚関係にある。

ヒント2（アプローチ）

制約をランダム順に並び替え、1本ずつ現在の最適点に追加する。満たすなら何もしない。満たさないなら新しい最適点は必ずその制約の境界線上にあるので、その直線上に限定した1次元LP（これまでの制約を直線上に制限したもの）を解けばよい。

ヒント3（誘導）

def line_1d_lp(con, active_cons, obj):
    a, b, c = con
    if abs(b) > 1e-12:
        P0, D = (0.0, c / b), (1.0, -a / b)
    else:
        P0, D = (c / a, 0.0), (0.0, 1.0)
    tlo, thi = -INF, INF
    for (aa, bb, cc) in active_cons:
        coef = aa * D[0] + bb * D[1]
        rhs = cc - (aa * P0[0] + bb * P0[1])
        if abs(coef) < 1e-12:
            if rhs < -1e-9: return None
        elif coef > 0:
            thi = min(thi, rhs / coef)
        else:
            tlo = max(tlo, rhs / coef)
    # tlo <= t <= thi の範囲でp*x+q*yを最小化するtを選ぶ

模範解答 (Python)

import sys
import random

INF = float('inf')


def solve_2d_lp(constraints, obj, box=1e7, seed=1):
    rnd = random.Random(seed)
    cs = constraints[:]
    rnd.shuffle(cs)

    p, q = obj
    box_cs = [(1, 0, box), (-1, 0, box), (0, 1, box), (0, -1, box)]

    vx = -box if p >= 0 else box
    vy = -box if q >= 0 else box
    v = (vx, vy)
    active = list(box_cs)

    def satisfies(pt, con):
        a, b, c = con
        return a * pt[0] + b * pt[1] <= c + 1e-9

    def line_1d_lp(con, active_cons):
        a, b, c = con
        if abs(b) > 1e-12:
            P0, D = (0.0, c / b), (1.0, -a / b)
        else:
            P0, D = (c / a, 0.0), (0.0, 1.0)
        tlo, thi = -INF, INF
        for (aa, bb, cc) in active_cons:
            coef = aa * D[0] + bb * D[1]
            rhs = cc - (aa * P0[0] + bb * P0[1])
            if abs(coef) < 1e-12:
                if rhs < -1e-9:
                    return None
            elif coef > 0:
                thi = min(thi, rhs / coef)
            else:
                tlo = max(tlo, rhs / coef)
        if tlo > thi + 1e-9:
            return None
        dcoef = p * D[0] + q * D[1]
        if dcoef > 0:
            t = tlo
        elif dcoef < 0:
            t = thi
        else:
            t = tlo if tlo > -INF else (thi if thi < INF else 0.0)
        if t in (INF, -INF):
            return None
        return (P0[0] + t * D[0], P0[1] + t * D[1])

    for con in cs:
        if satisfies(v, con):
            active.append(con)
            continue
        newv = line_1d_lp(con, active)
        if newv is None:
            return ('infeasible', None, None)
        v = newv
        active.append(con)

    if abs(abs(v[0]) - box) < 1e-6 or abs(abs(v[1]) - box) < 1e-6:
        return ('unbounded', None, None)
    val = p * v[0] + q * v[1]
    return ('optimal', v, val)


def solve():
    data = sys.stdin.read().split()
    idx = 0
    n = int(data[idx]); idx += 1
    p = float(data[idx]); q = float(data[idx + 1]); idx += 2
    constraints = []
    for _ in range(n):
        a = float(data[idx]); b = float(data[idx + 1]); c = float(data[idx + 2])
        idx += 3
        constraints.append((a, b, c))

    status, pt, val = solve_2d_lp(constraints, (p, q))
    if status == 'infeasible':
        print("INFEASIBLE")
    elif status == 'unbounded':
        print("UNBOUNDED")
    else:
        print(f"{val:.6f}")


solve()

Step-by-Step 解説

1仮想の箱制約で有界な問題に変換
最初から十分大きい箱|x|,|y|≤Mを制約として組み込み、常に有界な凸多角形の中でのLPとして扱う。初期最適点は箱の頂点として目的関数係数の符号だけで決まる。

2ランダム順序で制約を1本ずつ追加
制約の追加順序をランダムにすることで期待計算量の解析が成立する（バックワード解析）。

3違反したら1次元LPへ次元を落とす
新しい制約h_iを現在の最適点が満たさない場合、次の最適点は必ずh_iの境界線上にある。直線をP0+tDとパラメータ化し、これまでの制約をtに関する線形不等式に変換して[tlo,thi]の共通範囲で目的関数を最小化する。

4非有界判定
最終的な最適点が仮想の箱の境界に乗っていれば、本来の制約だけでは目的関数をいくらでも改善できることを意味しUNBOUNDEDと判定する。

5計算量
制約1本あたり満たされていればO(1)、違反していれば1次元LP。乱択増分法の解析により違反回数は期待O(log N)〜O(1)程度に抑えられ、全体で期待O(N)。

よくあるミス

ミス	原因	正しい書き方
制約をシャッフルせず入力順のまま処理する	乱択の意味を軽視し最悪ケースで計算量が悪化する	必ずrandom.shuffleなどで順序をランダム化してから処理する
直線のパラメータ化でb=0(垂直線)を特別扱いし忘れる	D=(1,-a/b)の式がゼロ除算を起こす	bがほぼ0のときはx=c/aでy方向にパラメータ化するなど場合分けを入れる
非有界判定を目的関数値の大きさで判定しようとする	箱のサイズと目的関数値の関係を混同する	最適点の座標が箱の境界に一致しているかで判定する
satisfiesの許容誤差を設けない	厳密な<=比較で浮動小数点の丸めにより誤判定する	ax+by <= c + eps のように小さな許容誤差を入れる

次のステップ

発展: 3次元以上のLP-type問題への一般化。次元dのとき乱択増分法は期待O(d!・N)程度で動く（dが定数なら線形）
次回予告: 次回のテーマは今後のセッションで選定

自己評価

完全理解おおむね理解要復習

Day 131-Q5 — Seidelのランダム化増分法（2次元線形計画法期待O(N)）

問題