F.Sondik — Alpha-vector value iteration in POMDPs

Notebook ID: 01.F.Sondik Area: 01 Decision-Making Under Uncertainty Status: Foundational Historical Thesis

Historical thesis under test: Smallwood & Sondik (1973) proved that the finite-horizon POMDP value function V_h(b) is piecewise linear and convex (PWLC) and can be represented exactly as the upper envelope of a finite set of linear functions over the belief simplex — the alpha-vectors. Each alpha-vector α ∈ ℝ^|S| defines a hyperplane α·b, and V_h(b) = max_α α·b. Bellman backup preserves this structure, so V_h is PWLC for all h if V_0 is.

Why this matters: The alpha-vector representation converts an uncountably infinite belief-state MDP into a finite (though growing) set of linear functions. This makes exact POMDP planning tractable for small state/action/observation spaces, and led directly to approximate methods (PBVI, SARSOP, HSVI) that scale to larger problems. It is the POMDP counterpart of Bellman's insight: optimal plans can be represented compactly, not enumerated.

Portfolio connection: F.Bayes showed how belief updates work (the filter). F.Shannon quantified how much each observation reduces uncertainty. F.Bellman showed that planning under uncertainty decomposes recursively. F.Sondik shows what the resulting value function looks like — and how to represent and compute it exactly. Together these four notebooks establish the foundation for everything in Area 01.

Setup

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
from scipy.stats import binom
from itertools import product

np.random.seed(42)

# -- Color palette (consistent with F-series) --
COL_TRUE   = "#2166AC"   # deep blue
COL_K1     = "#4393C3"
COL_K2     = "#92C5DE"
COL_K3     = "#D1E5F0"
COL_WARM   = "#D6604D"
COL_K5     = "#F4A582"
COL_GREY   = "#666666"
COL_GREEN  = "#1B7837"
COL_GOLD   = "#B8860B"

Tiger POMDP (same environment as F.Bellman)

State s ∈ {0=TigerLeft, 1=TigerRight} Action a ∈ {0=OpenLeft, 1=OpenRight, 2=Listen} Obs o ∈ {0=HearLeft, 1=HearRight} Belief b = P(tiger=Left) ∈ [0, 1]

Reward matrix R(s, a):

	s=Left	s=Right
OpenLeft	−100	+10
OpenRight	+10	−100
Listen	−1	−1

Transition: Opening resets to uniform (episode restart); Listen keeps state fixed. Observation model: P(HearLeft | TigerLeft, Listen) = p_listen = 0.85.

class TigerPOMDP:
    """Tiger POMDP (Cassandra et al. 1994 / Smallwood-Sondik testbed).

    All matrix representations needed for alpha-vector Bellman backup.
    """
    OPEN_L = 0;  OPEN_R = 1;  LISTEN = 2
    S = 2;       A = 3;       O = 2

    def __init__(self, p_listen: float = 0.85, gamma: float = 0.95):
        self.p  = p_listen
        self.gamma = gamma

        # R[s, a]  (reward received in state s taking action a)
        self.R = np.array([
            [-100.0,  10.0, -1.0],   # s=0 (TigerLeft)
            [  10.0,-100.0, -1.0],   # s=1 (TigerRight)
        ])

        # T[a, s, s']  (transition probabilities)
        # Listen: stays in same state.  Open: resets to uniform.
        self.T = np.zeros((self.A, self.S, self.S))
        for a in [self.OPEN_L, self.OPEN_R]:
            self.T[a, :, :] = 0.5   # reset to uniform
        self.T[self.LISTEN, 0, 0] = 1.0   # tiger stays
        self.T[self.LISTEN, 1, 1] = 1.0

        # Z[a, s', o]  (observation probabilities)
        # Only Listen is informative; open gives uniform obs (state is reset)
        self.Z = np.zeros((self.A, self.S, self.O))
        for a in [self.OPEN_L, self.OPEN_R]:
            self.Z[a, :, :] = 0.5
        p = self.p
        self.Z[self.LISTEN, 0, 0] = p;      self.Z[self.LISTEN, 0, 1] = 1.0 - p
        self.Z[self.LISTEN, 1, 0] = 1.0-p;  self.Z[self.LISTEN, 1, 1] = p

    def immediate_alpha(self, a: int) -> np.ndarray:
        """α₀^a[s] = R(s, a) — the base alpha-vector for action a."""
        return self.R[:, a].copy()

    def belief_from_b(self, b: float) -> np.ndarray:
        """Scalar b (= P(tiger=Left)) → probability vector [b, 1−b]."""
        return np.array([b, 1.0 - b])

    def dot_v(self, alphas: list, b_vec: np.ndarray) -> float:
        """V(b) = max_α  α · b  over the current set of alpha-vectors."""
        return max(np.dot(alpha, b_vec) for alpha in alphas)

    def best_alpha(self, alphas: list, b_vec: np.ndarray) -> np.ndarray:
        """Return the alpha-vector that achieves the maximum at b."""
        return max(alphas, key=lambda alpha: np.dot(alpha, b_vec))

Alpha-vector Bellman backup

Given the current set Γ_{h-1} of alpha-vectors, the Bellman backup produces a new set Γ_h. For each action a and each combination of one alpha-vector per observation (|O| choices), we compute one candidate alpha-vector:

$$\alpha^{a,\vec{k}}[s] = R(s,a) + \gamma \sum_{o} \sum_{s'} T(s'|s,a)\,Z(o|s',a)\,\alpha_{k_o}[s']$$

The result Γ_h is the subset of candidates that are ever optimal somewhere on the belief simplex (pruning removes dominated vectors).

For |S|=2, the belief simplex is the unit interval [0,1] and each alpha-vector is a line: α·b = α[0]·b + α[1]·(1−b). The upper envelope is easy to visualise.

def alpha_backup(pomdp: TigerPOMDP, prev_alphas: list) -> list:
    """One-step Bellman backup of an alpha-vector set.

    Returns the new (unpruned) set of candidate alpha-vectors.
    Each candidate corresponds to one (action, obs→alpha) combination.

    For the Tiger POMDP (|S|=2, |A|=3, |O|=2):
      - For each action a (3 choices)
      - For each pair (k0, k1) of alpha indices, one per obs (|Γ|² pairs)
    Total candidates before pruning: 3 × |Γ|² (grows exponentially without pruning).
    """
    g  = pomdp.gamma
    S  = pomdp.S
    A  = pomdp.A
    O  = pomdp.O
    candidates = []

    for a in range(A):
        # Precompute: for each (o, alpha_idx), the contribution vector
        # contrib[o][k][s] = Σ_{s'} T(s'|s,a)·Z(o|s',a)·α_k[s']
        contrib = np.zeros((O, len(prev_alphas), S))
        for o in range(O):
            for k, alpha_k in enumerate(prev_alphas):
                for s in range(S):
                    contrib[o, k, s] = sum(
                        pomdp.T[a, s, sp] * pomdp.Z[a, sp, o] * alpha_k[sp]
                        for sp in range(S)
                    )

        # For each tuple of alpha indices (k_o for each o)
        for idx_combo in product(range(len(prev_alphas)), repeat=O):
            alpha_new = pomdp.R[:, a].copy()   # immediate reward
            for o, k in enumerate(idx_combo):
                alpha_new += g * contrib[o, k]
            candidates.append(alpha_new)

    return candidates


def prune_dominated(alphas: list, n_points: int = 400) -> list:
    """Point-based pruning: keep only alpha-vectors that are best somewhere.

    For |S|=2 this is exact over a fine grid.  For larger state spaces,
    LP-based pruning (Lark's algorithm) is needed.
    """
    b_grid = np.linspace(0.0, 1.0, n_points)
    kept   = set()
    for i, b in enumerate(b_grid):
        bv  = np.array([b, 1.0 - b])
        best_val = -np.inf
        best_k   = 0
        for k, alpha in enumerate(alphas):
            v = np.dot(alpha, bv)
            if v > best_val:
                best_val = v
                best_k   = k
        kept.add(best_k)
    return [alphas[k] for k in sorted(kept)]


def value_iteration_alphas(pomdp: TigerPOMDP, max_h: int = 5,
                            n_prune: int = 400) -> list:
    """Compute alpha-vector sets Γ_0, Γ_1, ..., Γ_{max_h}.

    Γ_0: each action's immediate reward alpha-vector (3 vectors).
    Γ_h: Bellman backup of Γ_{h-1}, then prune.
    Returns list of sets indexed by horizon (0 to max_h).
    """
    # V_0(b) = 0 → represented as a single zero alpha-vector
    # (convention: start from zero value function)
    gamma_sets = [[ np.zeros(pomdp.S) ]]   # Γ_0 = {0}

    for h in range(1, max_h + 1):
        candidates = alpha_backup(pomdp, gamma_sets[-1])
        pruned     = prune_dominated(candidates, n_prune)
        gamma_sets.append(pruned)
        print(f"  h={h}: {len(candidates):4d} candidates → {len(pruned):3d} after pruning")

    return gamma_sets

E01 — Structure of the alpha-vector value function

Hypothesis: V_h(b) is piecewise linear and convex (PWLC) for all h. Each horizon adds alpha-vectors; the upper envelope grows in complexity but remains exactly representable by a finite number of linear pieces. At h=0, V_0(b) = 0 (one flat vector). At h=1 the three immediate reward lines dominate different regions. By h=3, the Listen action generates mixing effects that create intermediate pieces.

print("E01: Alpha-vector sets — PWLC structure")
print("=" * 50)

pomdp = TigerPOMDP(p_listen=0.85, gamma=0.95)
gamma_sets = value_iteration_alphas(pomdp, max_h=5)

b_fine = np.linspace(0.0, 1.0, 500)

horizons_plot = [1, 2, 3, 5]
fig, axes = plt.subplots(1, 4, figsize=(16, 4.5), sharey=False)
fig.suptitle("E01 — PWLC Value Function: Upper Envelope of Alpha-Vectors", y=1.01)

alpha_colors = [COL_TRUE, COL_K1, COL_WARM, COL_K5, COL_GREEN, COL_GOLD,
                COL_GREY, "#8B008B", "#2E8B57", "#FF8C00", "#00CED1", "#DC143C"]

for ax, h in zip(axes, horizons_plot):
    alphas = gamma_sets[h]
    V_env  = np.array([pomdp.dot_v(alphas, np.array([b, 1-b])) for b in b_fine])

    for k, alpha in enumerate(alphas):
        c = alpha_colors[k % len(alpha_colors)]
        line_vals = alpha[0] * b_fine + alpha[1] * (1.0 - b_fine)
        ax.plot(b_fine, line_vals, '--', color=c, lw=1.0, alpha=0.5)

    ax.plot(b_fine, V_env, '-', color='black', lw=2.5, label='V_h(b) envelope')

    ax.set_xlabel("b = P(tiger=Left)")
    ax.set_ylabel("Value")
    ax.set_title(f"h={h}  |Γ|={len(alphas)}")
    ax.set_xlim(0, 1)

    # annotate symmetry
    mid_v = np.interp(0.5, b_fine, V_env)
    ax.axvline(0.5, color='gray', ls=':', lw=1)
    ax.text(0.5, mid_v, f' V({0.5:.1f})={mid_v:.1f}', fontsize=7,
            color='gray', va='bottom')

plt.tight_layout()
plt.savefig('e01_sondik_pwlc.png', bbox_inches='tight')
plt.close()
print("\nSaved: e01_sondik_pwlc.png")

# Report alpha-vector counts
print("\nAlpha-vector set sizes:")
for h, gs in enumerate(gamma_sets):
    mid_v = pomdp.dot_v(gs, np.array([0.5, 0.5]))
    print(f"  h={h}: |Γ|={len(gs):3d}   V(0.5)={mid_v:.3f}")

FINDING 1: V_h(b) is PWLC at every horizon, confirmed by the upper-envelope structure. At h=1, three alpha-vectors define three linear segments, with the Listen vector dominant in the uncertain region near b=0.5 and the Open vectors dominant near certainty. Each successive backup adds vectors that capture finer distinctions between belief states. The value function is symmetric: V_h(b) = V_h(1−b) because the Tiger POMDP is symmetric in the two states. The alpha-vector set grows with h but remains finite — this is the Smallwood-Sondik theorem in action: finitely many vectors represent an uncountably infinite V_h.

E02 — Growth rate of alpha-vector sets

Hypothesis: The number of alpha-vectors |Γ_h| can grow exponentially with h in the worst case (one new vector per action-obs combination). In practice, pruning keeps the set small for structured problems. The Tiger POMDP's symmetry and the dominance of the Listen action's alpha-vectors in the uncertain region should keep |Γ_h| from exploding.

print("\nE02: Alpha-vector set growth")
print("=" * 50)

# Count vectors at each h and track pre/post pruning
b_test = np.array([0.5, 0.5])
results_e02 = {}

print(f"{'h':>4}  {'|Γ_h|':>8}  {'V(0.5)':>10}  {'V(0.1)':>10}  {'V(0.9)':>10}")
for h, gs in enumerate(gamma_sets):
    n_alpha = len(gs)
    v05 = pomdp.dot_v(gs, np.array([0.5, 0.5]))
    v01 = pomdp.dot_v(gs, np.array([0.1, 0.9]))
    v09 = pomdp.dot_v(gs, np.array([0.9, 0.1]))
    results_e02[h] = {'n': n_alpha, 'v05': v05, 'v01': v01, 'v09': v09}
    print(f"{h:>4}  {n_alpha:>8}  {v05:>10.3f}  {v01:>10.3f}  {v09:>10.3f}")

fig, axes = plt.subplots(1, 2, figsize=(12, 4.5))
fig.suptitle("E02 — Alpha-Vector Growth and Value Convergence", y=1.01)

ax = axes[0]
h_vals = list(results_e02.keys())
n_vals = [results_e02[h]['n'] for h in h_vals]
ax.bar(h_vals, n_vals, color=COL_TRUE, alpha=0.8)
for h, n in zip(h_vals, n_vals):
    ax.text(h, n + 0.1, str(n), ha='center', va='bottom', fontsize=10, fontweight='bold')
ax.set_xlabel("Horizon h")
ax.set_ylabel("|Γ_h| (number of alpha-vectors)")
ax.set_title("Alpha-vector set size vs horizon\n(Pruning prevents exponential blowup)")
ax.set_xticks(h_vals)

ax = axes[1]
for belief, col, lbl in [(0.1, COL_WARM, 'b=0.1 (near certainty left)'),
                          (0.5, COL_TRUE, 'b=0.5 (maximum uncertainty)'),
                          (0.9, COL_K1,  'b=0.9 (near certainty right)')]:
    key = {0.1: 'v01', 0.5: 'v05', 0.9: 'v09'}[belief]
    vals = [results_e02[h][key] for h in h_vals]
    ax.plot(h_vals, vals, 'o-', color=col, lw=2, ms=8, label=lbl)
ax.set_xlabel("Horizon h")
ax.set_ylabel("V_h(b)")
ax.set_title("Value convergence at three beliefs\n(Decreasing: listen cost dominates short horizons)")
ax.legend(fontsize=8)
ax.set_xticks(h_vals)

plt.tight_layout()
plt.savefig('e02_sondik_growth.png', bbox_inches='tight')
plt.close()
print("\nSaved: e02_sondik_growth.png")

FINDING 2: The alpha-vector set grows from 1 vector at h=0 to 13 at h=5, far below the theoretical maximum (3 × |Γ|² candidates per backup). Pruning is highly effective because the Tiger POMDP's symmetry and the Listen action's dominance in the central uncertain region eliminate most candidates. The set size is NOT monotone in h (h=4 has 7 vectors, fewer than h=3's 9) — this is correct: at some horizons the best action near the boundary can be achieved by fewer distinct vector segments.

V(0.5) is also non-monotone: dips to −1.95 at h=2 then jumps to +2.31 at h=3. This reflects a key threshold: h=2 cannot yet listen twice-then-open profitably (not enough steps), while h=3 CAN: the first profitable listen-listen-open sequence emerges exactly at h=3 for this environment. Near certainty (b=0.1, b=0.9), value is positive from h=2 onward because the agent can open immediately.

E03 — Alpha-vector identity: who "wins" where?

Hypothesis: Each alpha-vector corresponds to a specific policy fragment — a commitment to one action now and one alpha-vector from the previous horizon for each possible observation. The belief simplex is partitioned into regions where each alpha-vector (policy fragment) is optimal. Examining these regions reveals the policy structure: open-left dominates near b≈0, listen dominates near b=0.5, open-right dominates near b≈1.

print("\nE03: Alpha-vector regions at h=3 and h=5")
print("=" * 50)

b_fine2 = np.linspace(0.0, 1.0, 1000)

fig, axes = plt.subplots(2, 2, figsize=(14, 9))
fig.suptitle("E03 — Belief-Space Partition by Dominant Alpha-Vector", y=1.01)

for row, h in enumerate([3, 5]):
    alphas = gamma_sets[h]
    ax_val = axes[row, 0]
    ax_id  = axes[row, 1]

    V_vals  = []
    best_id = []

    for b in b_fine2:
        bv   = np.array([b, 1.0 - b])
        vals = [np.dot(alpha, bv) for alpha in alphas]
        V_vals.append(max(vals))
        best_id.append(int(np.argmax(vals)))

    # Value function
    ax_val.plot(b_fine2, V_vals, 'k-', lw=2.5)
    ax_val.set_xlabel("b = P(tiger=Left)")
    ax_val.set_ylabel("V_h(b)")
    ax_val.set_title(f"h={h}: Value function (PWLC)")

    # Color each alpha-vector
    colors_used = {}
    for k, alpha in enumerate(alphas):
        line = alpha[0]*b_fine2 + alpha[1]*(1.0-b_fine2)
        c = alpha_colors[k % len(alpha_colors)]
        ax_val.plot(b_fine2, line, '--', color=c, lw=1.0, alpha=0.4)
        colors_used[k] = c

    # Policy partition
    bid = np.array(best_id)
    ax_id.set_facecolor('#F5F5F5')
    for k in sorted(set(bid)):
        mask   = bid == k
        region = b_fine2[mask]
        if len(region) > 0:
            c = alpha_colors[k % len(alpha_colors)]
            ax_id.fill_between(region, 0, 1, color=c, alpha=0.6,
                               label=f"α_{k}")
            # Find midpoint of region
            mid = region[len(region)//2]
            ax_id.text(mid, 0.5, f'α_{k}', ha='center', va='center',
                       fontsize=9, fontweight='bold', color='white')

    ax_id.set_xlabel("b = P(tiger=Left)")
    ax_id.set_ylabel("")
    ax_id.set_title(f"h={h}: Dominant alpha-vector by belief region")
    ax_id.set_ylim(0, 1)
    ax_id.set_yticks([])
    ax_id.set_xlim(0, 1)

    # Mark boundaries
    boundaries = []
    for i in range(len(bid)-1):
        if bid[i] != bid[i+1]:
            boundaries.append(b_fine2[i])
            ax_id.axvline(b_fine2[i], color='black', lw=1.5, ls='-')
            ax_id.text(b_fine2[i], 0.05, f'{b_fine2[i]:.3f}',
                       ha='center', fontsize=7, rotation=90,
                       va='bottom', color='black')

    print(f"\nh={h}: {len(alphas)} alpha-vectors, {len(boundaries)} partition boundaries:")
    print(f"  Boundaries at b ≈ {[f'{x:.3f}' for x in boundaries]}")

plt.tight_layout()
plt.savefig('e03_sondik_regions.png', bbox_inches='tight')
plt.close()
print("\nSaved: e03_sondik_regions.png")

FINDING 3: The belief simplex is cleanly partitioned by the dominant alpha-vectors. At every horizon, the structure is symmetric: near b=0 (confident Tiger-is-Left), one set of alpha-vectors corresponding to OpenLeft is dominant; near b=1 (confident Tiger-is-Right), OpenRight-associated vectors dominate; in the middle uncertain region, Listen-associated vectors dominate. The partition boundaries shift with horizon: deeper lookahead sets a lower confidence threshold for acting, because the agent knows it can gather more information even after a suboptimal decision — but in the Tiger POMDP this effect is small because opens reset state. The PWLC partition is the dual view of the policy: it tells the agent which plan to commit to without enumerating all belief states.

E04 — Point-based approximation (PBVI)

The exact Smallwood-Sondik backup is tractable for the Tiger POMDP but exponential in |S| for larger problems. Point-Based Value Iteration (Pineau et al. 2003) selects a finite set B of reachable belief points and only backs up those points — producing one alpha-vector per point instead of one per action-obs combination. The resulting Γ is a subset of the exact Γ but covers the reachable belief space well in practice.

Hypothesis: PBVI with a representative set of belief points approximates the exact PWLC value function closely, with error concentrated at beliefs far from the sampled points. As |B| grows, the approximation tightens monotonically.

print("\nE04: Point-based value iteration vs exact backup")
print("=" * 50)

def pbvi_backup(pomdp: TigerPOMDP, prev_alphas: list,
                belief_points: np.ndarray) -> list:
    """PBVI backup: for each belief b in B, find the best backed-up alpha.

    For each b ∈ B and each action a:
      Compute the alpha that, when added to the action-specific contribution,
      maximises b·alpha_new.
    Keep one alpha per belief point (the one that maximises the value at that point).
    """
    S = pomdp.S
    O = pomdp.O
    g = pomdp.gamma
    new_alphas = []

    for b_vec in belief_points:
        best_val   = -np.inf
        best_alpha = None

        for a in range(pomdp.A):
            # Build alpha_a,b: greedy choice of one alpha per obs
            alpha_ab = pomdp.R[:, a].copy()
            for o in range(O):
                # For obs o under action a, find the alpha in prev_alphas
                # that maximises b'(b,a,o) · alpha  (approximately b · alpha
                # weighted by P(o|b,a))
                # Exact: evaluate each prev alpha at post-observation belief
                p_o_b = sum(
                    b_vec[s] * pomdp.T[a, s, sp] * pomdp.Z[a, sp, o]
                    for s in range(S) for sp in range(S)
                )
                if p_o_b < 1e-10:
                    # No information from this obs; pick zero contribution
                    continue

                # Contribution vector for this obs
                best_contrib_val = -np.inf
                best_contrib     = np.zeros(S)
                for alpha_k in prev_alphas:
                    contrib = np.zeros(S)
                    for s in range(S):
                        contrib[s] = sum(
                            pomdp.T[a, s, sp] * pomdp.Z[a, sp, o] * alpha_k[sp]
                            for sp in range(S)
                        )
                    # Evaluate at current b_vec (approximate; exact uses b')
                    val = np.dot(b_vec, contrib)
                    if val > best_contrib_val:
                        best_contrib_val = val
                        best_contrib     = contrib
                alpha_ab += g * best_contrib

            val = np.dot(b_vec, alpha_ab)
            if val > best_val:
                best_val   = val
                best_alpha = alpha_ab.copy()

        if best_alpha is not None:
            new_alphas.append(best_alpha)

    # Deduplicate (within tolerance)
    unique = []
    for alpha in new_alphas:
        if not any(np.allclose(alpha, u, atol=1e-8) for u in unique):
            unique.append(alpha)
    return unique


def pbvi_value(alphas: list, b: float) -> float:
    bv = np.array([b, 1.0 - b])
    return max(np.dot(alpha, bv) for alpha in alphas)


# Generate belief-point sets of different sizes
B_SIZES = [3, 7, 15, 30]
H_EVAL  = 5  # horizon to evaluate
EXACT_ALPHAS = gamma_sets[H_EVAL]

results_e04 = {}

print(f"Exact backup: {len(EXACT_ALPHAS)} alpha-vectors at h={H_EVAL}")

for n_B in B_SIZES:
    # Uniform belief points
    B_points = np.array([[b, 1.0-b] for b in np.linspace(0.0, 1.0, n_B)])

    # Run PBVI for H_EVAL steps from zero
    pb_alphas = [np.zeros(pomdp.S)]
    for h in range(1, H_EVAL + 1):
        pb_alphas = pbvi_backup(pomdp, pb_alphas, B_points)

    # Measure error vs exact
    b_eval = np.linspace(0.0, 1.0, 300)
    exact_vals = np.array([pomdp.dot_v(EXACT_ALPHAS, np.array([b, 1-b])) for b in b_eval])
    pbvi_vals  = np.array([pbvi_value(pb_alphas, b) for b in b_eval])
    errors = exact_vals - pbvi_vals   # PBVI always ≤ exact (lower bound property)

    max_err  = float(np.max(np.abs(errors)))
    mean_err = float(np.mean(np.abs(errors)))
    results_e04[n_B] = (pb_alphas, pbvi_vals, errors, max_err, mean_err)
    print(f"  |B|={n_B:3d}: {len(pb_alphas):3d} alpha-vectors, "
          f"max |error|={max_err:.3f}, mean |error|={mean_err:.3f}")

# Plot
fig, axes = plt.subplots(1, 2, figsize=(13, 5))
fig.suptitle(f"E04 — PBVI Approximation vs Exact Value Function (h={H_EVAL})", y=1.01)

ax = axes[0]
b_eval = np.linspace(0.0, 1.0, 300)
exact_vals = np.array([pomdp.dot_v(EXACT_ALPHAS, np.array([b, 1-b])) for b in b_eval])
ax.plot(b_eval, exact_vals, 'k-', lw=3, label='Exact', zorder=5)
pb_colors = [COL_WARM, COL_K5, COL_K1, COL_TRUE]
for (n_B, col) in zip(B_SIZES, pb_colors):
    ax.plot(b_eval, results_e04[n_B][1], '--', color=col, lw=1.5,
            label=f'PBVI |B|={n_B}')
ax.set_xlabel("b = P(tiger=Left)")
ax.set_ylabel("V_h(b)")
ax.set_title("PBVI approximates exact value function\n(Approaches exact as |B| grows)")
ax.legend(fontsize=8)

ax = axes[1]
for (n_B, col) in zip(B_SIZES, pb_colors):
    errs = results_e04[n_B][2]
    ax.plot(b_eval, errs, '-', color=col, lw=1.5, label=f'|B|={n_B}')
ax.axhline(0, color='k', ls='--', lw=1)
ax.set_xlabel("b = P(tiger=Left)")
ax.set_ylabel("V_exact(b) − V_PBVI(b)  [error]")
ax.set_title("Approximation error decays as |B| grows\n(Error ≥ 0: PBVI is a lower bound)")
ax.legend(fontsize=8)

plt.tight_layout()
plt.savefig('e04_sondik_pbvi.png', bbox_inches='tight')
plt.close()
print("\nSaved: e04_sondik_pbvi.png")

FINDING 4: PBVI produces a valid lower bound on the exact value function (V_PBVI ≤ V_exact everywhere). Error is highest at beliefs far from the sampled points. With |B|=3 (only at b=0, 0.5, 1), the approximation is coarse in between. With |B|=15–30, the error is small across the full simplex. The PBVI set is always a subset of the exact set in terms of quality: the same Bellman structure holds, but the coverage is restricted to the reachable belief space. This is the core insight of point-based methods: if the agent will only ever encounter beliefs near B, the exact computation for the rest of the simplex is wasted.

E05 — Convergence of the value function with horizon

Hypothesis: As h → ∞, V_h converges to the infinite-horizon optimal value function V* (for γ < 1). The Bellman contraction guarantees that ||V_h − V*||∞ ≤ γ^h · ||V_0 − V*||∞, so convergence is geometric. For the Tiger POMDP, the policy becomes stable (same dominant alpha-vector regions) well before the value function has fully converged.

print("\nE05: Value function convergence with horizon")
print("=" * 50)

# Compute alpha-sets up to h=10
print("Computing alpha-vector sets to h=10...")
gamma_sets_long = value_iteration_alphas(pomdp, max_h=10)

b_eval = np.linspace(0.0, 1.0, 300)
V_by_h = {}
for h, gs in enumerate(gamma_sets_long):
    V_by_h[h] = np.array([pomdp.dot_v(gs, np.array([b, 1-b])) for b in b_eval])

# Measure sup-norm difference between consecutive horizons
diffs = []
for h in range(1, len(gamma_sets_long)):
    d = float(np.max(np.abs(V_by_h[h] - V_by_h[h-1])))
    diffs.append(d)
    print(f"  ||V_{h} - V_{h-1}||_∞ = {d:.6f}")

# Geometric contraction rate
gamma_val = pomdp.gamma
print(f"\n  Discount factor γ = {gamma_val}")
print(f"  Theoretical contraction per step: γ = {gamma_val:.2f}")
if diffs:
    residual = diffs[-1]
    error_bound = gamma_val / (1 - gamma_val) * residual
    print(f"  Residual at h={len(diffs)}: ||V_{len(diffs)} − V_{len(diffs)-1}||_∞ = {residual:.4f}")
    print(f"  Convergence bound: γ/(1-γ) × residual = {gamma_val:.2f}/{1-gamma_val:.2f} × {residual:.4f} = {error_bound:.2f}")
    print(f"  → ||V_{len(diffs)} − V*||_∞ ≤ {error_bound:.2f}  (cannot claim V_h ≈ V*)")

fig, axes = plt.subplots(1, 2, figsize=(13, 5))
fig.suptitle("E05 — Convergence of V_h(b) to V*(b) (γ=0.95)", y=1.01)

ax = axes[0]
h_plot = [0, 1, 2, 3, 5, 7, 10]
cmap   = plt.cm.Blues(np.linspace(0.2, 1.0, len(h_plot)))
for h, c in zip(h_plot, cmap):
    if h < len(gamma_sets_long):
        ax.plot(b_eval, V_by_h[h], '-', color=c, lw=1.8, label=f'h={h}')
ax.set_xlabel("b = P(tiger=Left)")
ax.set_ylabel("V_h(b)")
ax.set_title("V_h(b) by horizon (h=0 to h=10; darkest = h=10)\n(Bellman contraction: consecutive gap → 0, but V_10 ≠ V*)")
ax.legend(fontsize=8, loc='upper right')

ax = axes[1]
h_range = list(range(1, len(gamma_sets_long)))
ax.semilogy(h_range, diffs, 'o-', color=COL_WARM, lw=2, ms=8, label='||V_h − V_{h-1}||_∞')
# Overlay theoretical rate
d0 = diffs[0]
theory = [d0 * (gamma_val ** (h-1)) for h in h_range]
ax.semilogy(h_range, theory, '--', color=COL_GREY, lw=1.5, label=f'γ^h · d₀  (γ={gamma_val})')
ax.set_xlabel("Horizon h")
ax.set_ylabel("||V_h − V_{h-1}||_∞  (log scale)")
ax.set_title("Geometric convergence of Bellman iterates\n(Actual ≤ theoretical bound)")
ax.legend(fontsize=9)
ax.set_xticks(h_range)

plt.tight_layout()
plt.savefig('e05_sondik_convergence.png', bbox_inches='tight')
plt.close()
print("\nSaved: e05_sondik_convergence.png")

FINDING 5: The Bellman iterates contract geometrically in sup-norm: each ratio ||V_{h+1} − V_h||∞ / ||V_h − V{h-1}||_∞ ≤ γ = 0.95, confirmed at every step. The ratio approaches 0.95 at h=4→5, showing the contraction bound is tight.

WHAT WE CANNOT CONCLUDE: that V_10 ≈ V*. The Bellman convergence theorem (Puterman 1994) gives: ||V_h − V*||∞ ≤ γ/(1-γ) · ||V_h − V{h-1}||∞. With γ=0.95 and residual ≈ 1.04 at h=10, the bound is 0.95/0.05 × 1.04 ≈ 19.8. A residual near 1 does NOT imply V* has been reached — the factor γ/(1-γ)=19 amplifies remaining error. Claiming V_10 ≈ V* requires demonstrating that the true gap ||V_10 − V*||∞ is small, not that consecutive iterates are close.

Importantly, pointwise values V_h(b) are NON-MONOTONE for problems with mixed-sign rewards. At b=0.5: V_0=0, V_1=−1, V_2=−1.95, V_3=+2.31, V_4=+1.80, V_5=+2.76. The value dips negative at h=1,2 (listen costs dominate), jumps positive at h=3 (now enough steps to listen twice and open profitably), then fluctuates. This is mathematically correct: the Bellman contraction guarantees convergence in sup-norm, not monotone convergence at each point.

Policy convergence may outpace value convergence — the partition boundaries from E03 show similar structure at h=3, 5, 7. However, we cannot confirm "essentially the infinite-horizon optimal policy" at h=5 without comparing against a known V* obtained by running many more iterations or by exact Lagrangian methods.

Summary

F.Sondik — Empirical conclusions:

| Experiment | Finding | |------------|---------| | E01 | V_h(b) is PWLC at every h; confirmed as upper envelope of finite alpha-vector set | | E02 | |Γ_h| non-monotone (pruning-dependent); V_h(0.5) oscillates before converging — h=3 is first horizon where listen-listen-open is profitable | | E03 | Belief simplex cleanly partitioned: Open-Left ≈ b<τ_L, Listen ≈ τ_L<b<τ_R, Open-Right ≈ b>τ_R | | E04 | PBVI is a valid lower bound; error ≥ 0 everywhere; tightens as |B| grows | | E05 | Geometric convergence in sup-norm (bound γ=0.95 tight at h=4); pointwise values non-monotone; policy stabilises before value function |

Theoretical result confirmed: Smallwood & Sondik (1973) — The finite-horizon POMDP value function is PWLC and exactly representable by a finite set of alpha-vectors. The upper-envelope structure V_h(b) = max_α α·b is preserved by every Bellman backup.

Connection to F-series sequence:

• F.Bayes: belief update (the filter that moves b_t → b_{t+1})
• F.Shannon: quantified information gain per observation
• F.Bellman: recursive decomposition of planning (Bellman equation for POMDPs)
• F.Sondik: representation of the resulting value function (PWLC, alpha-vectors)

These four together give a complete mathematical account of what it means for a system to know, decide, and act under partial observability — the foundational question of Area 01.

print("\nNotebook F.Sondik complete.")