The value of a state is an average you can solve for

One episode on a slippery grid

The grid in Fig. 1 has sixteen cells, numbered 0 to 15 row by row from the top left.

The environment takes an action and returns a reward and a next state. The agent is whatever supplies the action, and everything else is the environment.¹ A state is what the agent has in hand when it chooses, here its cell number. An action is one of up, right, down or left. A reward is the number the environment pays for one move. At time t the agent in state S_t takes action A_t, and one step later it receives reward R_{t+1} and is in state S_{t+1}. In Fig. 1 the agent in state 13 takes up, receives −1 and lands in state 9.

Moves slip. With slip probability 0.2, a move goes as intended with probability 0.8 and to each perpendicular side with probability 0.1, and a move into an edge leaves the agent where it is. Arriving at cell 3, the goal, pays +10. Arriving at cell 6, the pit, pays −20. Any other move pays −1. The goal and the pit are terminal states: arriving at one ends the run.¹ An episode is one run from a start state to a terminal state, seven moves in Fig. 1.

A policy, which Bellman's 1954 address framed as a prescription in terms of the current state, is the rule that picks the action in each state.² Formally it maps each state to a probability for each action, \pi(a \mid s).¹ The small arrows in Fig. 1 are one policy, drawn by hand, with probability 1 on a single action in each cell. From here on it stays fixed, and the question is how good each cell is under it.

The Markov decision process

The formal object is the Markov decision process, a decision process Bellman set out in 1957 and now written as the tuple (\mathcal{S}, \mathcal{A}, P, R, \gamma): states, actions, transition probabilities, rewards and a discount rate, the weight that shrinks later rewards.³⁴

For cell 5 and action up, P puts 0.8 on cell 1, 0.1 on cell 4 and 0.1 on cell 6. The reward here depends only on the arrival cell, so it is written R(s').

The Markov property is the assumption hidden in that notation: the next state depends on the history only through the current state and action.⁴

Read loosely, this says the past does not matter. The past does matter: it decided which cell the agent stands in. The property says that once the cell is known, the path that led there changes nothing about what comes next. It is a restriction on the state, which has to carry every aspect of the past that makes a difference to the future.¹

Value is an average of returns

The return G_t is the sum of the rewards after time t, each weighted by a power of the discount rate \gamma.¹

Fig. 1 uses \gamma = 0.9. That episode returned 0.63: six moves at −1 with falling weights, then the goal's +10 at weight 0.53. The same sum obeys a recursion, G_t = R_{t+1} + \gamma G_{t+1}.¹

Any five-move walk to the goal, the fewest moves possible, returns 3.12, and 8,386 of the episodes did that. The sample mean is −1.94 with a standard error of 0.05. "How good is cell 13" is taken to mean this average, because slips make episodes from one cell differ and no single episode answers it.

The value of state s under policy \pi is that average taken over all episodes, the expected return from s, and it is 0 at a terminal state.¹ Exactly, v_\pi(13) = -2.01.

The Bellman expectation equation

The definition needs whole episodes. Peeling off the first move gives an equation that needs only one move. The first step replaces the return with its recursion.

The next step averages over the arrival cell. Let P_\pi(s' \mid s) = \sum_a \pi(a \mid s) P(s' \mid s, a) be the probability that the policy and the environment together take s to s'. Landing in s' fixes the first reward at R(s').

This is the step that needs the Markov property. Given S_{t+1} = s', each later state depends only on the state and action before it, and the policy picks each action from the current state alone. So everything later has a distribution that does not depend on S_t, and the condition on S_t drops. Neither the transition probabilities nor the policy changes with time, so what is left is the definition of v_\pi(s').

Substituting gives the Bellman expectation equation.⁴

The intended move contributes 2.107 and the 0.1 chance of the pit contributes −2.000, so a cell three moves from the goal is worth 0.06. The arithmetic uses v_\pi(1), and v_\pi(1) has v_\pi(5) in its own equation, so one substitution does not finish the job.

Cycles between cells force a joint solve

The Bellman expectation equation is one linear equation for each of the 14 non-terminal cells, in 14 unknowns.¹ Fig. 4 shows which unknowns each equation uses.

At slip 0 each cell has one successor and the graph has no cycles, so the equations solve back to front in one pass: v(2) = 10, then v(1) = -1 + 0.9 \times 10 = 8, and on through the order in Fig. 4.⁵ The code runs that pass, and it matches a matrix solve to within 9 \times 10^{-16}.

Slip 0.2 adds arrows that rule this out. Cell 1 can slip down to cell 5, and cell 5 moves up to cell 1, so each value appears in the other's equation and neither can be computed first. A cell appearing in its own equation is a smaller matter. Cell 2 can slip up into the edge and stay put, which only takes a division.

Only cycles through two different cells force a joint solve. If the cells can be ordered so that each equation uses only itself and earlier cells, then at \gamma < 1 one pass with one division per cell solves the system. Such an order exists exactly when no cycle runs through two or more cells, because two cells on one cycle would each have to come before the other.

Randomness is not the cause: turning cell 1 to face cell 0 at slip 0 makes a cycle with no randomness in it, and cells 0 and 1 then have to be solved together. At slip 0.2 the groups in Fig. 4 are the graph's strongly connected components, and solving them in order, each from the ones before it, also matches the matrix solve.⁶

The direct solve

Let Q be the matrix P_\pi restricted to non-terminal cells, since the goal and the pit have value 0. Let b(s) = \sum_{s'} P_\pi(s' \mid s) R(s') be the expected reward of one move from s. All the equations together read v = b + \gamma Q v, which solves in one step.⁴

Sweeps

Iterative policy evaluation starts from a guess and applies the equation as an update, v_{k+1} = b + \gamma Q v_k, replacing every cell once per sweep.¹

From v_0 = 0, sweep k holds the expected discounted reward over the first k moves from each cell, so news of the goal and the pit spreads outward one cell per sweep. Cell 12 is six moves from the goal. It sinks to −4.43 at sweep 5 while only step costs and the pit are in reach, and has climbed to −2.19 by sweep 8. The converged panel shows the risk: cell 13 is one move closer to the goal than cell 12 and is worth less, −2.01 against −1.86, because its episodes spend more moves beside the pit. The straight route from 13 passes the pit at cells 5 and 2, and a slip to the right takes it past cells 10 and 7. The straight route from 12 passes the pit only at cell 2.

Why sweeps converge under a discount

Write T(v) = b + \gamma Q v and \lVert x \rVert_\infty = \max_s |x(s)|, the max norm. For \gamma < 1 the map T is a contraction in this norm: it multiplies the distance between any two vectors by at most \gamma.⁷⁸

The inequality holds because each entry of Qx is a sum of entries of x with non-negative weights totalling at most 1. The same bound gives \lVert (\gamma Q)^k \rVert_\infty \le \gamma^k, so \sum_k (\gamma Q)^k converges. Multiplying by I - \gamma Q telescopes to I, so the sum is the inverse, and v exists and is unique. Since T(v) = v, putting u = v_k and w = v gives one factor of \gamma per sweep, from every start.

Only a policy that never ends needs the discount

The discount sets how much a reward counts by how far ahead it lies: k moves ahead, it carries weight \gamma^k. At γ = 0 each cell shows its expected next reward. Cell 0, three moves from the goal, climbs through −1.00, −0.52, 2.16 and 3.43 as γ goes from 0 to 1.

The last panel has no discount at all, and the solve and the sweeps both work: the largest error is 1.35 \times 10^{-7} after 30 sweeps. The contraction proof gives nothing here, because its factor is 1 and cell 0's row of Q sums to 1. What works instead is that this policy is proper: every cell has a path of positive-probability moves to a terminal cell.⁹ Take m long enough for all those paths and let \epsilon be the smallest probability, over cells, of ending within m moves, which is positive because each cell has such a path. Row sums of Q^m are survival probabilities, so \lVert Q^m \rVert_\infty \le 1 - \epsilon and \lVert Q^{jm} \rVert_\infty \le (1-\epsilon)^j \to 0.¹⁰ Since \lVert Q \rVert_\infty \le 1, the powers in between are no larger, so Q^k \to 0. That puts the spectral radius \rho(Q), the largest eigenvalue magnitude, below 1 and makes \sum_k Q^k the inverse of I - Q, so v exists and is unique.¹⁰ The error after k sweeps is Q^k (v_0 - v), so it goes to 0 from every start.

Changing four arrows so that the bottom row presses against the bottom wall gives the policy called STUCK in the code. A slip there only moves along the wall, so no episode that enters the bottom row ever leaves it.

The bottom row's block of Q has rows summing to 1, so \rho(Q) = 1 and I - Q has no inverse. Every move there pays −1, so at γ = 1 cell 13 falls by 1 per sweep without end. Every γ below 1 restores a finite answer: paying −1 forever is worth -1/(1-\gamma), which is −10 at γ = 0.9 and −100 at γ = 0.99.¹ With that −1 per move, a cell falls exactly when an episode from it can get stuck, so at γ = 1 cells 8 to 15 fall and the top two rows settle at the values in Fig. 6. A discount gives every policy an answer, and a terminal cell gives one only to policies that reach it. Whether cells must be solved together is a separate question: the fixed policy at slip 0.2 has cycles and still ends with probability 1.

The measured rate

The error after k sweeps is (\gamma Q)^k (v_0 - v), and the asymptotic factor per sweep of such an iteration is the spectral radius of its matrix.¹¹ Here \rho(\gamma Q) = 0.4404, measured at 0.4406 by sweep 30, against the 0.9 of the bound, which is the worst case and which the stuck policy attains. The radius has to be that of Q, restricted to non-terminal cells, because the goal and the pit hold value 0 and carry no error. Kept as absorbing rows, they add an eigenvalue 1 for every policy, which would put the radius of the scaled matrix at exactly γ. Its eigenvectors are nonzero at the goal or the pit, so sweeps from v_0 = 0 never excite it. The guide in Fig. 8 is anchored at sweep 30 and describes the tail only, and the first sweeps fall more slowly than it.

Two files that print the results

With both files in one folder, uv run evaluate.py installs numpy from the inline metadata and finishes in about a second. Changing the line GAMMA, SLIP = 0.9, 0.2 reruns everything at another setting. gridworld.py holds the environment and the two policies as arrays. The sweep keeps two arrays, the old values and the new ones.¹

# /// script
# requires-python = ">=3.10"
# dependencies = ["numpy"]
# ///
"""A slippery 4x4 grid as arrays. Cells are 0 to 15, row by row from top."""

import numpy as np

N = 16
GOAL, PIT = 3, 6
TERMINAL = [GOAL, PIT]
ARROWS = "^>v<"  # the actions up, right, down, left are numbered 0, 1, 2, 3
STEPS = [(-1, 0), (0, 1), (1, 0), (0, -1)]

# R[s'] is the reward for arriving in s'.
R = np.full(N, -1.0)
R[GOAL], R[PIT] = 10.0, -20.0


def draw(rows):
    """A policy drawn as rows of arrows. Terminal cells, ".", get action 0."""
    return np.array([max(ARROWS.find(arrow), 0) for arrow in rows.split()])


# The fixed policy, written by hand.
PI = draw("""> > > .
             ^ ^ . ^
             ^ ^ > ^
             ^ ^ ^ ^""")

# The same policy with the bottom row pressed against the bottom wall.
STUCK = draw("""> > > .
                ^ ^ . ^
                ^ ^ > ^
                v v v v""")


def move(s, a):
    """The cell reached by moving as intended. An edge leaves s in place."""
    row, col = s // 4 + STEPS[a][0], s % 4 + STEPS[a][1]
    return 4 * row + col if 0 <= row < 4 and 0 <= col < 4 else s


def transitions(slip):
    """P[s, a, s'], with all-zero rows for the terminal cells."""
    P = np.zeros((N, 4, N))
    for s in set(range(N)) - set(TERMINAL):
        for a in range(4):
            P[s, a, move(s, a)] += 1 - slip
            P[s, a, move(s, (a + 1) % 4)] += slip / 2
            P[s, a, move(s, (a + 3) % 4)] += slip / 2
    return P


def fixed_policy(policy, slip):
    """Q, the policy's moves among non-terminal cells, and b, the expected
    one-step reward from each cell. Terminal rows and columns of Q are zero."""
    P_pi = transitions(slip)[np.arange(N), policy]
    Q = P_pi.copy()
    Q[:, TERMINAL] = 0.0
    return Q, P_pi @ R


def episode(P, policy, start, rng):
    """One episode as a list of (cell, intended action, reward, next cell)."""
    s, steps = start, []
    while s not in TERMINAL:
        a = policy[s]
        s_next = rng.choice(N, p=P[s, a])
        steps.append((s, a, R[s_next], s_next))
        s = s_next
    return steps


def discounted(steps, gamma):
    """The return of a list of steps, from its first step."""
    return sum(gamma**t * step[2] for t, step in enumerate(steps))


def radius(A):
    """The spectral radius: the largest eigenvalue magnitude."""
    return abs(np.linalg.eigvals(A)).max()

# /// script
# requires-python = ">=3.10"
# dependencies = ["numpy"]
# ///
"""Evaluate fixed policies on the grid and print the page's results."""

import numpy as np
from gridworld import (N, PI, PIT, R, STUCK, TERMINAL, discounted, episode,
                       fixed_policy, radius, transitions)

GAMMA, SLIP = 0.9, 0.2
START, EPISODES = 13, 20_000
LIVE = [s for s in range(N) if s not in TERMINAL]

def solve(Q, b, gamma):
    return np.linalg.solve(np.eye(N) - gamma * Q, b)

def sweeps(Q, b, gamma, k):
    vs = [np.zeros(N)]
    for _ in range(k):
        vs.append(b + gamma * Q @ vs[-1])
    return np.array(vs)

def components(Q):
    """Groups of cells that all reach each other, each after those it enters"""
    reach = np.linalg.matrix_power(np.eye(N) + (Q > 0), N) > 0
    groups = {tuple(t for t in LIVE if reach[s, t] and reach[t, s])
              for s in LIVE}
    return sorted(groups, key=lambda g: (reach[g[0]].sum(), g))

def solve_by_components(Q, b, gamma):
    v = np.zeros(N)
    for g in map(list, components(Q)):
        v[g] = np.linalg.solve(np.eye(len(g)) - gamma * Q[np.ix_(g, g)],
                               b[g] + gamma * Q[g] @ v)
    return v

def show(title, *grids):
    print(title)
    for row in range(4):
        print("   ".join(" ".join(f"{x:5.2f}" for x in g[4 * row:][:4])
                         for g in grids))

def main():
    P, (Q, b) = transitions(SLIP), fixed_policy(PI, SLIP)
    v = solve(Q, b, GAMMA)
    steps = episode(P, PI, START, np.random.default_rng(2))
    print("episode", START, *[s[3] for s in steps], "return so far",
          *np.cumsum([GAMMA**t * s[2] for t, s in enumerate(steps)]).round(2))
    rng = np.random.default_rng(0)
    runs = [episode(P, PI, START, rng) for _ in range(EPISODES)]
    G = np.array([discounted(run, GAMMA) for run in runs])
    pit = np.mean([run[-1][3] == PIT for run in runs])
    print(f"mean {G.mean():.3f} +- {G.std() / EPISODES**0.5:.3f}, exact "
          f"{v[START]:.3f}, pit {pit:.3f}, goal {1 - pit:.3f}, "
          f"top {G.max():.2f}")
    print(f"unit bins from -20:\n{np.histogram(G, np.arange(-20, 7))[0]}")
    branch = P[5, PI[5]] * (R + GAMMA * v)
    print("cell 5:", *branch[branch != 0].round(3), f"sum {branch.sum():.3f}, "
          f"v(1) {v[1]:.3f}, v(4) {v[4]:.3f}")
    turned = np.where(np.arange(N) == 1, 3, PI)  # cell 1 faces left, to cell 0
    for name, policy, slip in (("slip 0", PI, 0.0), (f"slip {SLIP}", PI, SLIP),
                               ("slip 0, cell 1 turned", turned, 0.0)):
        Qs, bs = fixed_policy(policy, slip)
        # 0.9, not GAMMA: cell 1 turned never ends, so gamma 1 has no solution.
        gap = abs(solve_by_components(Qs, bs, 0.9) - solve(Qs, bs, 0.9))
        print(f"{name}, gap {gap.max():.0e}:",
              *[g if g[1:] else g[0] for g in components(Qs)])
    vs = sweeps(Q, b, GAMMA, 30)
    show("sweeps 1, 2, 3", vs[1], vs[2], vs[3])
    show("sweeps 5, 8, exact", vs[5], vs[8], v)
    show("gamma 0, 0.5, 1", *[solve(Q, b, g) for g in (0.0, 0.5, 1.0)])
    for g in (GAMMA, 1.0):
        e = abs(sweeps(Q, b, g, 30) - solve(Q, b, g)).max(axis=1)
        print(f"gamma {g}: rho {radius(g * Q):.4f}, sweep 30 ratio "
              f"{e[30] / max(e[29], 1e-99):.4f}, error {e[30]:.2e}, "
              f"bound {g**30 * e[0]:.2e}")
    full = P[np.arange(N), PI] + np.diag(np.isin(np.arange(N), TERMINAL))
    print(f"absorbing rows kept, gamma {GAMMA}: "
          f"radius {radius(GAMMA * full):.4f}")
    Qn, bn = fixed_policy(STUCK, SLIP)
    vs = sweeps(Qn, bn, 1.0, 4000)
    finite = [round(solve(Qn, bn, g)[13]) for g in (0.9, 0.99)]
    print(f"stuck: rho(Q) {radius(Qn):.4f}, v(13), k = 1000, gamma 0.9, 0.99:",
          round(vs[1000, 13]), *finite)
    print("falling at 4000:", *np.flatnonzero(abs(vs[4000] - vs[3999]) > 1e-9))

if __name__ == "__main__":
    main()

episode 13 9 5 4 0 1 2 3 return so far -1.0 -1.9 -2.71 -3.44 -4.1 -4.69 0.63
mean -1.944 +- 0.053, exact -2.005, pit 0.197, goal 0.803, top 3.12
unit bins from -20:
[   0 1605    0 1246  708  287   88   12    0    0    0    0    0    4
   11   68  302  336  883 1223 3138 1703    0 8386    0    0]
cell 5: 2.107 -0.045 -2.0 sum 0.062, v(1) 4.037, v(4) 0.613
slip 0, gap 9e-16: 2 7 1 11 0 5 10 15 4 9 14 8 13 12
slip 0.2, gap 9e-16: 2 7 (0, 1, 4, 5) (8, 9, 10, 11, 12, 13, 14, 15)
slip 0, cell 1 turned, gap 4e-15: 2 7 (0, 1) 11 4 5 10 15 8 9 14 12 13
sweeps 1, 2, 3
-1.00 -1.00  5.90  0.00   -1.90  2.90  6.43  0.00    0.73  3.56  6.48  0.00
-1.00 -2.90  0.00  5.90   -2.07 -3.71  0.00  6.43   -2.89 -1.00  0.00  6.48
-1.00 -1.00 -2.90 -1.00   -1.90 -3.44 -3.71  2.90   -2.97 -4.18 -1.11  3.56
-1.00 -1.00 -1.00 -1.00   -1.90 -1.90 -3.27 -1.90   -2.71 -3.94 -4.01  0.62
sweeps 5, 8, exact
 1.85  3.96  6.48  0.00    2.13  4.03  6.48  0.00    2.16  4.04  6.48  0.00
-0.14 -0.17  0.00  6.48    0.55  0.04  0.00  6.48    0.61  0.06  0.00  6.48
-2.12 -1.83 -0.29  3.95   -0.86 -1.10 -0.13  4.02   -0.72 -1.03 -0.10  4.02
-4.43 -3.03 -1.81  1.72   -2.19 -2.20 -1.14  1.95   -1.86 -2.01 -1.08  1.98
gamma 0, 0.5, 1
-1.00 -1.00  5.90  0.00   -0.52  1.44  6.21  0.00    3.43  4.85  6.56  0.00
-1.00 -2.90  0.00  5.90   -1.40 -2.40  0.00  6.21    2.07  1.18  0.00  6.56
-1.00 -1.00 -2.90 -1.00   -1.76 -2.17 -2.43  1.43    0.74  0.12  0.96  4.82
-1.00 -1.00 -1.00 -1.00   -1.90 -2.07 -2.10 -0.56   -0.56 -0.96 -0.01  3.17
gamma 0.9: rho 0.4404, sweep 30 ratio 0.4406, error 5.24e-09, bound 2.75e-01
gamma 1.0: rho 0.4894, sweep 30 ratio 0.4895, error 1.35e-07, bound 6.56e+00
absorbing rows kept, gamma 0.9: radius 0.9000
stuck: rho(Q) 1.0000, v(13), k = 1000, gamma 0.9, 0.99: -1000 -10 -100
falling at 4000: 8 9 10 11 12 13 14 15