Day 079-Q4 — 確率的最短路（SSP・ベルマン方程式 + ガウス消去法）

2026-07-01 赤色 Master / Phase 8+ ★★★★★★★★★ SSP・マルコフ連鎖・期待値・ガウス消去・線形方程式系

問題

$N$ 頂点 $M$ 辺の有向グラフがある。各頂点 $i$ にはコスト $c_i \ge 0$ が付いており、確率的遷移を行う。頂点 $i$ から確率 $p_{ij}$ で頂点 $j$ へ移動しコスト $c_i$ を支払う。ゴール頂点 $t$ に到達するまでのコストの期待値 $E[s]$ を求めよ。

$$E[i] = c_i + \sum_{j=1}^{N} p_{ij} \cdot E[j] \quad (i \neq t), \quad E[t] = 0$$

制約

パラメータ範囲備考
$N$$2 \le N \le 200$頂点数
$M$$1 \le M \le N^2$辺数
$c_i$$0 \le c_i \le 10^6$各頂点のコスト
$p_{ij}$有理数（分子/分母で入力）各頂点からの確率和=1

入出力例

入力例1

4 4 1 4
5 3 2 0
1 2 1 2
1 3 1 2
2 4 1 1
3 4 1 1

出力例1

15/2

$E[4]=0$, $E[2]=3$, $E[3]=2$, $E[1]=5+\frac{1}{2}(3)+\frac{1}{2}(2)=\frac{15}{2}$

概念図: 確率グラフ → ベルマン方程式 → ガウス消去

確率的最短路: ベルマン方程式 → 連立一次方程式 → ガウス消去 確率グラフ（入力例1） s=1 c=5 2 c=3 3 c=2 t=4 c=0 p=1/2 p=1/2 p=1 p=1 E[1] = ? E[2]=3, E[3]=2 E[4]=0 → E[1] = 5 + 1/2·3 + 1/2·2 = 15/2 連立方程式 → ガウス消去 ベルマン方程式（E[4]=0 代入後）: E[1] − 1/2·E[2] − 1/2·E[3] = 5 1·E[2] = 3 1·E[3] = 2 拡大係数行列 [A | b]: │ 1 −1/2 −1/2 │ 5 │ │ 0 1 0 │ 3 │ │ 0 0 1 │ 2 │ ↓ ガウス・ジョルダン消去 │ 1 0 0 │15/2│ │ 0 1 0 │ 3 │ │ 0 0 1 │ 2 │ → E[s] = E[1] = 15/2

各頂点で「コスト + 次頂点の期待値の加重平均」というベルマン方程式を立て、移項して線形方程式系に変換。 ガウス消去法で厳密解（有理数）を求める。

ヒント

ヒント1（方向性）

ゴール頂点 $t$ 以外の各頂点 $i$ に対して未知数 $E[i]$ を立てる。$E[i] - \sum_{j \neq t} p_{ij} E[j] = c_i$ という線形方程式系 $AE = b$ をガウス消去法で $O(N^3)$ で解く。

ヒント2（アプローチ）

係数行列 $A$（$(N-1) \times (N-1)$）の設定:
対角 $A[i][i] = 1$、非対角 $A[i][j] = -p_{v_i, v_j}$（遷移確率）、定数 $b[i] = c_{v_i}$。
fractions.Fraction を使うと厳密な有理数解が得られる。

ヒント3（ほぼ答え） 
from fractions import Fraction
# mat: (n x n+1) 拡大係数行列
for col in range(n):
    pivot = next((r for r in range(col, n) if mat[r][col] != 0), -1)
    if pivot == -1: continue
    mat[col], mat[pivot] = mat[pivot], mat[col]
    inv = Fraction(1, mat[col][col])
    for k in range(n + 1):
        mat[col][k] *= inv
    for row in range(n):
        if row != col and mat[row][col] != 0:
            f = mat[row][col]
            for k in range(n + 1):
                mat[row][k] -= f * mat[col][k]
# 解: mat[i][n] = E[v_i]

模範解答

import sys
from fractions import Fraction

def solve():
    data = sys.stdin.read().split()
    ptr = 0
    N = int(data[ptr]); ptr += 1
    M = int(data[ptr]); ptr += 1
    s = int(data[ptr]); ptr += 1
    t = int(data[ptr]); ptr += 1

    cost = [0] * (N + 1)
    for v in range(1, N + 1):
        cost[v] = int(data[ptr]); ptr += 1

    edges_from = [[] for _ in range(N + 1)]
    for _ in range(M):
        u  = int(data[ptr]); ptr += 1
        v  = int(data[ptr]); ptr += 1
        pn = int(data[ptr]); ptr += 1
        pd = int(data[ptr]); ptr += 1
        edges_from[u].append((v, Fraction(pn, pd)))

    nodes   = [v for v in range(1, N + 1) if v != t]
    n       = len(nodes)
    idx_map = {v: i for i, v in enumerate(nodes)}

    mat = [[Fraction(0)] * (n + 1) for _ in range(n)]
    for i, v in enumerate(nodes):
        mat[i][i] = Fraction(1)
        mat[i][n] = Fraction(cost[v])
        for (w, p) in edges_from[v]:
            if w != t:
                j = idx_map[w]
                mat[i][j] -= p

    for col in range(n):
        pivot_row = next((r for r in range(col, n) if mat[r][col] != 0), -1)
        if pivot_row == -1: continue
        mat[col], mat[pivot_row] = mat[pivot_row], mat[col]
        inv_p = Fraction(1, mat[col][col])
        for k in range(n + 1):
            mat[col][k] *= inv_p
        mat[col][col] = Fraction(1)
        for row in range(n):
            if row != col:
                f = mat[row][col]
                if f != 0:
                    for k in range(col, n + 1):
                        mat[row][k] -= f * mat[col][k]
                    mat[row][col] = Fraction(0)

    E = {t: Fraction(0)}
    for i, v in enumerate(nodes):
        E[v] = mat[i][n]

    ans = E[s]
    p, q = ans.numerator, ans.denominator
    print(p if q == 1 else f"{p}/{q}")

solve()

Step-by-Step 解説

Step 1: ベルマン方程式の導出

頂点 $i$ での期待総コスト $E[i]$ を定義する。ゴール: $E[t] = 0$。その他: $$E[i] = c_i + \sum_j p_{ij} E[j]$$ 直感: まずコスト $c_i$ を支払い、その後 $j$ に移動したときの残り期待コスト $E[j]$ の加重平均を加える。 サイクルがあると後ろ向きDPが使えないため、連立方程式として解く。

Step 2: 線形方程式系への変換

$E[t]=0$ を代入して移項: $E[i] - \sum_{j \neq t} p_{ij} E[j] = c_i$

$A$ の要素値理由
$A[i][i]$$1$$E[i]$ の係数
$A[i][j]$（$j \neq i, j \neq t$）$-p_{v_i, v_j}$$-p_{ij} \cdot E[j]$ の係数
$b[i]$$c_{v_i}$定数項（右辺）

行列 $A$ は対角優位（$|A_{ii}|=1 \ge \sum_{j \neq i}|A_{ij}|$）なので数値的に安定。

Step 3: ガウス・ジョルダン消去

拡大係数行列 $[A | b]$ に対して各列をピボットとして前進・後退消去を同時実行（ジョルダン型）。終了時に $[I | E_{\text{vec}}]$ の形になり直接解が得られる。計算量 $O(N^3)$。$N=200$ なら $8 \times 10^6$ 演算で十分高速。

Step 4: fractions.Fraction による厳密解

浮動小数点では累積誤差が生じる。fractions.Fraction を使うと有理数を厳密に扱え、最終的な期待値を真の既約分数として出力できる。 $N=200$ の場合 Fraction の演算量は許容範囲内。出力: $Q=1$ なら整数、そうでなければ P/Q 形式。

Step 5: 計算量まとめ

ステップ計算量
方程式設定$O(NM)$
ガウス消去$O(N^3)$
全体$O(N^3) = O(200^3) = 8 \times 10^6$

よくあるミス

ミス原因正しい書き方
ゴール頂点 $t$ も変数に含める 方程式が $(N \times N)$ になりランク不足 $t$ を変数から除き $E[t]=0$ を代入してから立式
対角 $A[i][i]=1$ の初期化漏れ $E[i]$ の係数が0になり方程式が $-\sum p_{ij}E[j]=c_i$ になる 必ず mat[i][i] = 1 を設定
float を使う 分数の累積誤差で整数解が 7.4999... になる fractions.Fraction を使う
ゴールへの遷移 $p_{it}$ を定数項に加算 $-p_{it} \cdot E[t] = 0$ なので不要だが加算するとバグ ゴールへの遷移は何もしない
ジョルダン消去せず前進消去のみ 後退代入のコードが別途必要で複雑 全行を対象に消去するジョルダン型が実装簡潔
到達不能頂点を見落とす ピボットが0になった列でゼロ除算が起きる if pivot_row == -1: continue でスキップ

次のステップ

  • 発展問題1: コストが負になり得る場合（MDPの Value Iteration と収束証明）
  • 発展問題2: 複数ゴールが存在し、最速で任意ゴールに到達するまでの期待ステップ数
  • 発展問題3: 無限ループが存在し得るグラフで $E[s] = \infty$ を検出して報告する（ランク不足の判定）
  • 発展問題4: $N=1000$ 以上の大規模版で浮動小数点 Value Iteration を使う場合の収束判定

自己評価