Exact Solution Methods

Chapter 1: Policy Evaluation — What Is a Policy Worth?

Before we can optimize, we need to evaluate. Given a fixed policy π, how good is it? The answer is the value function U^π(s): the expected discounted return starting from s and following π forever.

Computing U^π is called the prediction problem. It has a clean recursive structure. From any state s, following policy π gives immediate reward R(s, π(s)), then we land in state s′ with probability T(s′|s, π(s)), and get value U^π(s′) from there. Putting it together, the Bellman expectation equation (eq. 7.6) must hold at every state simultaneously:

This is a system of |S| equations in |S| unknowns. Two ways to solve it:

Why Does the Iterative Method Converge?

The update is a contraction mapping. Define the Bellman operator B^π as “apply one round of the update to every state.” For any two value functions U₁ and U₂:

After one application, the gap between any two guesses shrinks by at least a factor of γ. Because γ < 1, repeated application squeezes the gap to zero — regardless of starting point. The fixed point they converge to is exactly U^π.

Iterative Evaluation: Worked Trace

Consider the 3-state chain MDP from Chapter 0 (states 0, 1, 2; actions left/right). Policy π: always go right. γ = 0.9. R(1, right) = 10, all else −1.

Iteration 0: U = [0, 0, 0].

Iteration 1: State 0: R(0,right)=−1, T(1|0,right)=1, so U(0)=−1+0.9×0=−1. State 1: R(1,right)=10, T(2|1,right)=1, so U(1)=10+0.9×0=10. State 2: stays at 2 (wall), R=−1, U(2)=−1.

Iteration 2: U(0)=−1+0.9×10=8. U(1)=10+0.9×(−1)=9.1. U(2)=−1+0.9×(−1)=−1.9.

Iteration 3: U(0)=−1+0.9×9.1=7.19. U(1)=10+0.9×(−1.9)=8.29. U(2)=−1+0.9×(−1.9)=−2.71.

Values are propagating backward from the reward. State 0 is slowly learning it’s one step away from a big payoff. The exact fixed point (via matrix inversion) is U^π = [−3.45, 6.36, −5.91] — which is what the iteration converges to over many steps. (The negative values at states 2 reflect the “dead end” of hitting the wall repeatedly.)

python
import numpy as np

def iterative_policy_evaluation(mdp, pi, k_max=100, eps=1e-6):
    """Evaluate policy pi by repeated Bellman updates (Algorithm 7.3)."""
    U = {s: 0.0 for s in mdp.S}
    for _ in range(k_max):
        U_new = {}
        for s in mdp.S:
            a = pi[s]
            U_new[s] = mdp.R(s, a) + mdp.gamma * sum(
                mdp.T(s, a, sp) * U[sp] for sp in mdp.S
            )
        delta = max(abs(U_new[s] - U[s]) for s in mdp.S)
        U = U_new
        if delta < eps:
            break
    return U

def exact_policy_evaluation(mdp, pi):
    """Exact matrix solve: U = (I - gamma*T_pi)^{-1} R_pi (Algorithm 7.4)."""
    n = len(mdp.S)
    idx = {s: i for i, s in enumerate(mdp.S)}
    R_pi = np.array([mdp.R(s, pi[s]) for s in mdp.S])
    T_pi = np.array([[mdp.T(s, pi[s], sp) for sp in mdp.S] for s in mdp.S])
    U_vec = np.linalg.solve(np.eye(n) - mdp.gamma * T_pi, R_pi)
    return {s: U_vec[idx[s]] for s in mdp.S}

Chapter 2: Greedy Policies and Q-Functions

Suppose someone hands you a value function U. Maybe it’s an estimate, maybe it’s exact. Either way, you want to extract the best policy from it. The natural idea: at each state, pick the action that looks best over the next step.

Formally, the greedy policy with respect to U (eq. 7.11) is:

If U = U* (optimal), the greedy policy is optimal. If U is approximate, the greedy policy is approximately optimal — and we can bound exactly how suboptimal it is.

The Q-Function

The expression inside the argmax has a name: the action value function (or Q-function), defined in eq. 7.12:

Q(s,a) is the expected return of: (1) taking action a in state s, then (2) acting greedily forever after (using value function U to decide). From Q, we get everything for free:

The Advantage Function

The advantage function A(s,a) (eq. 7.15) measures how much better action a is compared to the best action:

Since U(s) = max_a Q(s,a), the greedy action has A(s, π(s)) = 0. Every other action has A(s,a) ≤ 0 — a negative advantage tells you exactly how much you lose by choosing that action.

Worked Example: Q-Table From U

Here is a toy 2×3 grid MDP with γ = 0.9. Goal state (bottom-right) has been solved; its value is U*(goal) = 0. All move costs are −1. Suppose the value function at state (0,0) is U(0,0) = −5, and from there you can move right to (0,1) with U(0,1) = −4, or down to (1,0) with U(1,0) = −4.

Q(s, right) = −1 + 0.9×(−4) = −1 − 3.6 = −4.6.

Q(s, down) = −1 + 0.9×(−4) = −4.6 (tie; both neighbors equally far from goal).

Q(s, left) = −1 + 0.9×U(wall=stay) = −1 + 0.9×(−5) = −5.5 (wall, loops to same state).

Q(s, up) = −1 + 0.9×U(wall=stay) = −5.5 (same; another wall).

So U(s) = max Q = −4.6. A(s, right) = −4.6−(−4.6) = 0. A(s, left) = −5.5−(−4.6) = −0.9. Going left wastes 0.9 units of expected return.

python
def lookahead(mdp, U, s, a):
    """Q(s,a): one-step lookahead value (Algorithm 7.2)."""
    return mdp.R(s, a) + mdp.gamma * sum(
        mdp.T(s, a, sp) * U[sp] for sp in mdp.S
    )

def greedy(mdp, U, s):
    """Extract greedy action at state s."""
    return max(mdp.A, key=lambda a: lookahead(mdp, U, s, a))

def Q_table(mdp, U):
    """Full Q-table: Q[s][a] for all states and actions."""
    return {s: {a: lookahead(mdp, U, s, a) for a in mdp.A} for s in mdp.S}

def advantage(mdp, U, s, a):
    """A(s,a) = Q(s,a) - U(s). Always <= 0."""
    return lookahead(mdp, U, s, a) - U[s]

Chapter 3: Policy Iteration

We know how to evaluate a policy (Chapter 1) and how to improve it (Chapter 2). The obvious question: what if we alternated? Evaluate, then improve, then evaluate again — until the policy stops changing?

That’s exactly policy iteration (Algorithm 7.6). Starting from any policy π, it cycles:

Why Does It Converge?

Three facts lock together to guarantee convergence:

(1) Monotonic improvement. The policy improvement theorem says: if π′ is greedy with respect to U^π, then U^π′(s) ≥ U^π(s) for all states. The new policy is never worse. Proof: at each state, π′ picks the best one-step lookahead, which is at least as large as what π would have gotten. Chain this inequality forward through all future states (just like we did in Chapter 2) and you get U^π′(s) ≥ U^π(s) globally.

(2) Finiteness. A finite MDP with |S| states and |A| actions has at most |A|^|S| deterministic policies (each state independently chooses one of |A| actions). Each improvement step either finds a strictly better policy or terminates. No cycles are possible — we can’t visit the same policy twice, because each visit strictly improves it.

(3) Optimality at termination. If π′ = π, then the greedy policy equals the current one, which means U^π(s) = max_a[R(s,a) + γ∑T(s′|s,a)U^π(s′)] for every state. That IS the Bellman optimality equation. So U^π = U* and π = π*.

The Modified Policy Iteration Spectrum

Using full matrix evaluation at each step is expensive. What if we only run a few sweeps of iterative evaluation before improving? This gives modified policy iteration:

Hex World Example (Figure 7.5 from the book)

Start with the policy “always go east.” Iteration 1: evaluate (values propagate from reward state), improve (greedy extraction changes some arrows near the reward). Iteration 2: values update again, more arrows change. Iteration 3: converged. Only 3 full iterations for the exact optimal policy.

The early iterations change many arrows. Later iterations change only a few. Convergence is front-loaded — the policy improves fastest at the start, then makes tiny refinements near the end.

python
def policy_iteration(mdp, pi_init=None, k_eval=50):
    """Algorithm 7.6: Policy iteration."""
    # Initialize to first action if not provided
    pi = pi_init or {s: mdp.A[0] for s in mdp.S}
    while True:
        # 1. Evaluate current policy
        U = iterative_policy_evaluation(mdp, pi, k_max=k_eval)
        # 2. Improve greedily
        pi_new = {s: greedy(mdp, U, s) for s in mdp.S}
        # 3. Check convergence
        if all(pi_new[s] == pi[s] for s in mdp.S):
            return pi, U
        pi = pi_new

Chapter 4: Value Iteration (Showcase)

What if we skipped the separate evaluation step entirely? Instead of first computing U^π and then extracting the greedy policy, we combine them into a single update. At every state, we immediately apply the optimal Bellman equation — the max over all actions — and use that as the new value estimate.

This is the Bellman backup (eq. 7.16), the engine of value iteration:

Starting from U(s) = 0 everywhere and applying this update to every state repeatedly, we converge to the optimal value function U*, which satisfies the Bellman optimality equation (eq. 7.17):

Why Does It Converge?

The max operator does not break the contraction. For any two value functions U₁ and U₂, applying the Bellman optimality operator T*:

The key identity is |max_a f(a) − max_a g(a)| ≤ max_a |f(a) − g(a)|. Both players pick their best action; the gap between their best values can’t exceed the worst-case gap across all actions. And that worst-case gap is bounded by γ times the original difference in U₁, U₂. Since γ < 1, repeated application squeezes any initial error to zero.

Stopping and Error Bounds

We can’t measure ||U_k − U*|| directly (we don’t know U*). Instead, we track the Bellman residual:

The true error is bounded by: ||U_k − U*||_∞ ≤ δγ/(1−γ). And the policy loss — how much worse the extracted policy is than optimal — is bounded by: ||U^π − U*||_∞ ≤ 2εγ/(1−γ) where ε is the value error.

Worked Example: 3-State MDP by Hand

States: s₀, s₁, s₂ (chain). Actions: stay (s), advance (a). R(s₁, advance) = 10; all else −1. Advance from s₀ goes to s₁; advance from s₁ goes to s₂; s₂ is absorbing (loops). Discount γ = 0.9.

k=0: U = [0, 0, 0].

k=1:
 s₀: max(stay: −1+0.9×0=−1, advance: −1+0.9×0=−1) = −1
 s₁: max(stay: −1+0.9×0=−1, advance: 10+0.9×0=10) = 10
 s₂: max(stay: −1+0.9×0=−1, advance: −1+0.9×0=−1) = −1

k=2:
 s₀: max(−1+0.9×(−1)=−1.9, −1+0.9×10=8) = 8 ← advance is now clearly better!
 s₁: max(−1+0.9×10=8, 10+0.9×(−1)=9.1) = 9.1
 s₂: max(−1+0.9×(−1)=−1.9, −1+0.9×(−1)=−1.9) = −1.9

k=3:
 s₀: advance: −1+0.9×9.1=7.19
 s₁: advance: 10+0.9×(−1.9)=8.29
 s₂: — = −2.71

Already at k=2, value iteration “learned” that s₀ should advance (reaching s₁ gives 10). The backup propagates optimism: once s₁’s value is updated to 10, that knowledge immediately ripples to s₀ on the next step.

python
def bellman_backup(mdp, U, s):
    """One Bellman backup at state s (Algorithm 7.7)."""
    return max(
        mdp.R(s, a) + mdp.gamma * sum(mdp.T(s, a, sp) * U[sp] for sp in mdp.S)
        for a in mdp.A
    )

def value_iteration(mdp, eps=1e-6, k_max=1000):
    """Algorithm 7.8: Value iteration with Bellman residual stopping."""
    U = {s: 0.0 for s in mdp.S}
    for k in range(k_max):
        U_new = {s: bellman_backup(mdp, U, s) for s in mdp.S}
        delta = max(abs(U_new[s] - U[s]) for s in mdp.S)
        U = U_new
        if delta < eps * (1 - mdp.gamma) / mdp.gamma:
            # Tighter stop: guarantees true error < eps
            break
    policy = {s: greedy(mdp, U, s) for s in mdp.S}
    return policy, U

Chapter 5: Asynchronous Value Iteration

Standard value iteration is synchronous: in each iteration, we compute new values for ALL states simultaneously, using only the old values. Then we swap in the new table and repeat. Clean, correct — but it ignores useful information. When we compute U(s₅), we already have a fresh U(s₁). Why not use it?

Asynchronous value iteration (Algorithm 7.9) updates states one at a time, in place, using the most recent value of every state already updated this sweep:

The difference: when we update U(s₅), we use the already-fresh U(s₁), U(s₂), etc. This is the Gauss-Seidel method, adapted from linear algebra to dynamic programming.

State Ordering Matters

The order in which we update states controls how fast information flows. There are three natural orderings:

One more key fact: even with a bad ordering, asynchronous VI cannot converge slower than synchronous VI. It can only be the same or faster. The only requirement for convergence: every state must be updated infinitely often. Any ordering that visits every state repeatedly works.

Worked Comparison: Same Problem, Two Orderings

Linear chain of 5 states (s₀→s₁→s₂→s₃→s₄). Reward at s₄: +10. γ = 0.9. All transitions deterministic (always advance).

Left-to-right ordering (s₀, s₁, s₂, s₃, s₄): Sweep 1 updates s₀ first — it sees U(s₁)=0, stays at 0. Then s₁ sees U(s₂)=0 → stays 0. … Only s₄ has the reward. After 1 sweep: same as synchronous VI. Needs ~5 sweeps to fully propagate.

Right-to-left ordering (s₄, s₃, s₂, s₁, s₀): Sweep 1 updates s₄ first: U(s₄) = 10. Then s₃ sees fresh U(s₄)=10: U(s₃) = −1 + 0.9×10 = 8. Then s₂ sees fresh U(s₃)=8: U(s₂) = 6.2. s₁: 4.58. s₀: 3.12. Full propagation in 1 sweep!

python
def gauss_seidel_vi(mdp, ordering=None, k_max=1000, eps=1e-6):
    """Algorithm 7.9: Gauss-Seidel (asynchronous) value iteration."""
    # ordering: list of states in update order. Default: mdp.S as-is.
    ordering = ordering or mdp.S
    U = {s: 0.0 for s in mdp.S}
    for _ in range(k_max):
        delta = 0.0
        for s in ordering:           # update in-place, use fresh values immediately
            u_new = bellman_backup(mdp, U, s)
            delta = max(delta, abs(u_new - U[s]))
            U[s] = u_new              # write back immediately (in-place)
        if delta < eps:
            break
    return {s: greedy(mdp, U, s) for s in mdp.S}, U

# Backward induction: sweep from last time step to first
def backward_induction(mdp_h, H):
    """Finite-horizon DP: sweep backward from t=H. Returns policy for each t.
    mdp_h.S assumed to contain (state, time) tuples.
    """
    U = {(s, H): mdp_h.terminal_reward(s) for s in mdp_h.base_states}
    policies = {}
    for t in range(H - 1, -1, -1):    # sweep backward
        for s in mdp_h.base_states:
            best_a, best_v = None, -float('inf')
            for a in mdp_h.A:
                v = mdp_h.R(s, a) + mdp_h.gamma * sum(
                    mdp_h.T(s, a, sp) * U.get((sp, t+1), 0)
                    for sp in mdp_h.base_states
                )
                if v > best_v: best_a, best_v = a, v
            U[(s, t)] = best_v
            policies[(s, t)] = best_a
    return policies, U

Chapter 6: Linear Program Formulation

Value iteration and policy iteration are iterative algorithms — they converge asymptotically. Is there a way to compute U* in one shot, as a direct optimization? Yes: rewrite the Bellman optimality equation as a linear program.

The Bellman optimality equation requires U*(s) = max_a[R(s,a) + γ∑T(s′|s,a)U*(s′)] for every state s. The “max” is the hard part — it makes the equation nonlinear. The LP trick: replace the max with a constraint. Instead of saying “U(s) equals the best lookahead,” we say “U(s) is at least as large as every lookahead”:

This is a system of |S|×|A| linear constraints in |S| variables. The feasible set is a convex polyhedron. But which point in this polyhedron should we pick? Any U satisfying the constraints is an overestimate of U*. The smallest feasible U is exactly U*.

So the LP is (eq. 7.19–7.20):

Complexity: MDPs Are in P

The LP has |S| variables and |S|×|A| constraints. LP is solvable in polynomial time (via interior-point methods). Therefore, MDPs with finite state and action spaces are solvable in polynomial time — they lie in the complexity class P.

A 2-State Example by Hand

Two states s₀, s₁. Two actions: left (L), right (R). γ = 0.9. Transitions: L always stays; R from s₀ goes to s₁; R from s₁ stays. Rewards: R(s₀,L)=−1, R(s₀,R)=−1, R(s₁,L)=−1, R(s₁,R)=10.

LP constraints (variables: U₀, U₁):

From s₀, action L: U₀ ≥ −1 + 0.9 U₀ → 0.1 U₀ ≥ −1 → U₀ ≥ −10

From s₀, action R: U₀ ≥ −1 + 0.9 U₁

From s₁, action L: U₁ ≥ −1 + 0.9 U₁ → 0.1 U₁ ≥ −1 → U₁ ≥ −10

From s₁, action R: U₁ ≥ 10 + 0.9 U₁ → 0.1 U₁ ≥ 10 → U₁ ≥ 100

The binding constraint on U₁ is U₁ = 100. Then U₀ ≥ −1 + 0.9×100 = 89. Minimize: U₀ = 89, U₁ = 100. The optimal policy: in s₀, go right (worth 89); in s₁, go right (keeps collecting reward 10, discounted sum = 100).

python
from scipy.optimize import linprog
import numpy as np

def lp_solve(mdp):
    """Algorithm 7.10: Solve MDP via linear programming."""
    S, A, gamma = mdp.S, mdp.A, mdp.gamma
    n = len(S)
    idx = {s: i for i, s in enumerate(S)}

    # Objective: minimize sum of U(s)
    c = np.ones(n)

    # Constraints: U(s) >= R(s,a) + gamma * sum_s' T(s'|s,a) U(s')
    # Rearranged: -U(s) + gamma * sum_s' T(s'|s,a) U(s') <= -R(s,a)
    rows = []
    rhs = []
    for s in S:
        for a in A:
            row = np.zeros(n)
            row[idx[s]] = -1.0
            for sp in S:
                row[idx[sp]] += gamma * mdp.T(s, a, sp)
            rows.append(row)
            rhs.append(-mdp.R(s, a))

    res = linprog(c, A_ub=np.array(rows), b_ub=np.array(rhs), method='highs')
    U = {s: res.x[idx[s]] for s in S}
    return {s: greedy(mdp, U, s) for s in S}, U

Chapter 7: Linear Quadratic Regulator

Every method so far assumed a finite state and action space. But most real physical systems — a pendulum, a drone, a robotic arm — have continuous states and actions. You can’t enumerate them. Can we still compute an exact optimal policy?

Yes — for one important special case: when dynamics are linear and rewards are quadratic. This is the Linear Quadratic Regulator (LQR), and it has a closed-form solution.

The Setup

Linear dynamics (eq. 7.22): The next state is a linear function of the current state and action, plus Gaussian noise:

Here T_s is the state transition matrix (n×n), T_a is the action transition matrix (n×m), and w is zero-mean Gaussian noise.

Quadratic reward (eq. 7.23): The reward penalizes both large states and large actions:

R_s is negative semidefinite (penalizes large states); R_a is negative definite (penalizes large actions). The goal: keep the system near zero with minimal control effort.

The Solution

The optimal value function turns out to be quadratic: U_h(s) = s^TV_hs + q_h. The optimal policy turns out to be linear: π_h(s) = L_h · s.

The gain matrix L_h (eq. 7.29):

V_h satisfies the discrete-time Riccati equation (eq. 7.30):

Worked Example 7.4: Position-Velocity Control

State: s = [x, v]. Action: scalar acceleration a. Timestep Δt = 0.1. Dynamics: x′ = x + vΔt + ½aΔt², v′ = v + aΔt. Reward: R(s,a) = −x² − v² − 0.5a².

In matrix form: T_s = [[1, 0.1],[0, 1]], T_a = [[0.005],[0.1]], R_s = diag(−1,−1), R_a = [−0.5].

Starting at s₀ = [−10, 0], the Riccati recursion gives gains:

With a longer horizon the controller pushes harder early: at x=−10, h=5 gives initial action ≈ +5.04, rapidly accelerating toward the origin, while h=2 gives only +2.86.

python
import numpy as np

def lqr(Ts, Ta, Rs, Ra, h_max=50):
    """Algorithm 7.11: LQR via backward Riccati iteration."""
    n = Ts.shape[0]
    V = np.zeros((n, n))    # V_0 = 0 at terminal step
    gains = []
    for _ in range(h_max):
        M = Ta.T @ V @ Ta + Ra
        L = -np.linalg.solve(M, Ta.T @ V @ Ts)   # eq 7.29
        gains.append(L)
        V = Rs + Ts.T @ V @ Ts + Ts.T @ V @ Ta @ L  # Riccati, eq 7.30
    return gains   # gains[h-1] = L_h

# Example 7.4 setup
dt = 0.1
Ts = np.array([[1.0, dt], [0.0, 1.0]])
Ta = np.array([[0.5*dt**2], [dt]])
Rs = np.diag([-1.0, -1.0])
Ra = np.array([[-0.5]])
gains = lqr(Ts, Ta, Rs, Ra, h_max=10)
print("L_5 =", gains[4])    # [[-0.504, -1.124]]

# Simulate trajectory from s0 = [-10, 0]
s = np.array([-10.0, 0.0])
L = gains[-1]   # use converged gain
trajectory = [s.copy()]
for _ in range(50):
    a = L @ s        # linear policy: a = L*s
    s = Ts @ s + Ta.reshape(2) * a[0]
    trajectory.append(s.copy())

Chapter 8: Value Iteration Sandbox

You’ve seen value iteration on fixed examples. Now you’re in control. Below is a fully configurable grid world: choose the size, place rewards and walls by clicking, set the discount factor and stochastic slip probability. Run value iteration step by step or let it animate to convergence. Watch the Bellman residual decay. Extract and display optimal policy arrows.

Chapter 9: Policy vs. Value Iteration Side by Side

Both methods solve the same problem and both produce the same answer: the optimal policy. But they take very different paths to get there. The comparison below runs them simultaneously on an identical MDP so you can see the tradeoff directly.

Per-Iteration Cost Breakdown

For a grid with N = |S| states and A = |A| actions:

For the 5×5 grid above (N=25, A=4): PI costs 25³ ≈ 15,000 operations per outer iter but typically finishes in 3 outer iters. VI costs 25²×4 = 2,500 operations per iter but needs ∼50–200 iters for γ=0.9. Total: PI ≈ 45,000 vs. VI ≈ 125,000–500,000. PI wins for small grids.

Chapter 10: Method Comparison

Five exact methods, each with its own niche. Here’s the map:

When Each Method Breaks

The Curse of Dimensionality

All exact methods share a fatal flaw: they require enumerating all states. A robot with 10 joints, each discretized to 100 positions, has 100¹⁰ = 10²⁰ states. No computer can enumerate them. This is the curse of dimensionality — the fundamental barrier that drives the entire field of approximate dynamic programming (Chapter 8) and model-based RL (Chapter 9).

Chapter 11: Connections and What’s Next

Chapter 7 is the foundation. Everything that follows in the book is an answer to one of three problems that exact methods cannot handle: state spaces that are too large, models that are unknown, and partial observability.

Equation Cheat Sheet

Related Lessons

These micro-lessons complement Chapter 7:

• microLearning: Markov Decision Processes — visual intro to MDPs and value iteration from scratch
• microLearning: RL Algorithms — Q-learning, SARSA, policy gradient when the model is unknown

“An optimal policy has the property that whatever the initial state and initial decision are, the remaining decisions must constitute an optimal policy with regard to the state resulting from the first decision.”
— Richard Bellman, Dynamic Programming, 1957