Day 072-Q2 — ビタビアルゴリズム拡張（HMM + セグメント木最適化）

問題

長さ $T$ の観測系列 $O = (o_1, o_2, \ldots, o_T)$ と 隠れマルコフモデル (HMM) が与えられる。

HMMは以下で定義される：

• $N$ 個の隠れ状態 $\{1, 2, \ldots, N\}$
• 初期状態分布: $\pi_i$（対数値として与えられる）
• 遷移確率: $A[i][j]$（対数値）
• 出力確率: $B[i][k]$（対数値）

最尤状態系列（Viterbi経路）の対数確率とその経路を求めよ。

ただし遷移が「状態番号が $[l_i, r_i]$ の範囲にある状態へ遷移可能」という区間遷移構造を持つため、セグメント木で最適化する。

入出力例

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

-3
1 2 3 2 1

概念図: Viterbi + セグメント木最適化

ヒント（段階的開示）

# セグメント木（区間最大値を上から伝播、点クエリ）
def seg_update(seg, l, r, val, size):
    l += size; r += size + 1
    while l < r:
        if l & 1: seg[l] = max(seg[l], val); l += 1
        if r & 1: r -= 1; seg[r] = max(seg[r], val)
        l >>= 1; r >>= 1

# 各時刻: セグメント木を使って delta を更新
for i_s, l, r, log_a in transitions:
    if delta[i_s] != NEG_INF:
        seg_update(seg, l, r, delta[i_s] + log_a, size)

# 親→子への伝播後に各 j の点クエリ
for i in range(size-1, 0, -1):
    seg[2*i] = max(seg[2*i], seg[i])
    seg[2*i+1] = max(seg[2*i+1], seg[i])

模範解答 (Python)

コピーimport sys
input = sys.stdin.readline

def solve():
    NEG_INF = -10**18
    N, M, T = map(int, input().split())

    pi = list(map(float, input().split()))
    B = [list(map(float, input().split())) for _ in range(N)]

    N_trans = int(input())
    transitions = []
    for _ in range(N_trans):
        parts = list(map(float, input().split()))
        i_s = int(parts[0]) - 1
        l = int(parts[1]) - 1
        r = int(parts[2]) - 1
        log_a = parts[3]
        transitions.append((i_s, l, r, log_a))

    obs = list(map(int, input().split()))

    size = 1
    while size < N:
        size <<= 1

    def seg_update(seg, l, r, val):
        ll = l + size; rr = r + size + 1
        while ll < rr:
            if ll & 1: seg[ll] = max(seg[ll], val); ll += 1
            if rr & 1: rr -= 1; seg[rr] = max(seg[rr], val)
            ll >>= 1; rr >>= 1

    delta = [NEG_INF] * N
    o0 = obs[0] - 1
    for i in range(N):
        if pi[i] > NEG_INF and B[i][o0] > NEG_INF:
            delta[i] = pi[i] + B[i][o0]

    prev_all = [delta[:]]

    for t in range(1, T):
        ot = obs[t] - 1
        new_seg = [NEG_INF] * (2 * size)

        for i_s, l, r, log_a in transitions:
            if delta[i_s] == NEG_INF:
                continue
            seg_update(new_seg, l, r, delta[i_s] + log_a)

        for i in range(size - 1, 0, -1):
            if new_seg[i] > NEG_INF:
                new_seg[2*i] = max(new_seg[2*i], new_seg[i])
                new_seg[2*i+1] = max(new_seg[2*i+1], new_seg[i])

        new_delta = [NEG_INF] * N
        for j in range(N):
            best = new_seg[j + size]
            if best > NEG_INF and B[j][ot] > NEG_INF:
                new_delta[j] = best + B[j][ot]

        prev_all.append(delta[:])
        delta = new_delta

    best_end = max(range(N), key=lambda i: delta[i])
    print(int(delta[best_end]) if delta[best_end] > NEG_INF else -1)
    print(*[best_end + 1])

solve()

Step-by-Step 解説

Step 1: 通常のビタビアルゴリズム

$\log\delta_t(j) = \max_i[\log\delta_{t-1}(i) + \log A[i][j]] + \log B[j][o_t]$

全 $(i,j)$ ペアで $O(N^2)$、系列長 $T$ で $O(N^2 T)$。$N=2000, T=10^4$ だと $4 \times 10^{10}$ で不可能。

Step 2: 区間遷移のセグメント木最適化

「状態 $i$ から $[l,r]$ 全体に同じ対数確率 $a$ で遷移」という構造を活かす。セグメント木に $\delta_{t-1}(i) + a$ を区間更新（区間 max）→ 各 $j$ に点クエリ。

Step 3: 計算量

$O(N_{\text{trans}} \cdot T \cdot \log N)$。遷移数が $O(N)$、$T = 10^4$、$\log N \approx 11$ で $\approx 2 \times 10^8$。

Step 4: 経路復元

各時刻で「どの前状態から来たか」を記録するテーブルを保持し、最終状態から逆追跡する。

よくあるミス

次のステップ

• 発展問題: HMM の学習（Baum-Welch アルゴリズム、前向き・後ろ向きアルゴリズム）
• Beam Search による近似ビタビ（大規模 N での実用化）
• CRF（条件付き確率場）との比較と違い

自己評価