Decisions
Solving a finite Markov decision process given as tensors — value iteration, exact policy evaluation, policy iteration and advantages.
Prelude.Decisions solves a finite Markov decision process (MDP) when you know its dynamics. An MDP is two tensors: p of shape S×A×S, where p[s, a, t] is the probability that action \(a\) in state \(s\) leads to state \(t\), and r of shape S×A, the expected reward of taking \(a\) in \(s\). Everything else is a contraction against those two and a maximum over the action axis.
Terminal states are written as absorbing states with zero reward: p[s, a, s] = 1 and r[s, a] = 0 for every action. Every state keeps the same number of actions, so shapes stay fixed, and a terminal value is pinned at zero by construction.
The number of sweeps is an argument, not a tolerance. Hilbert has no data-dependent loop count, so you ask for \(k\) sweeps and read drift to see whether \(k\) was enough. The number needed to reach accuracy \(\varepsilon\) is about \(\ln(\varepsilon(1-\gamma)/(2 r_{\max})) / \ln \gamma\); 200 sweeps at \(\gamma = 0.9\) is far past machine precision.
Exact evaluation uses the normal equations. solve is conjugate gradient, which requires a symmetric positive-definite matrix, and \(I - \gamma P^\pi\) is not symmetric. The module therefore solves \(A^\top A\,v = A^\top r^\pi\), which is valid for every \(\gamma < 1\) but squares the condition number: \(\kappa(I - \gamma P^\pi) \le (1+\gamma)/(1-\gamma)\), which is 19 at \(\gamma = 0.9\) and 199 at 0.99. Beyond \(\gamma = 0.999\) prefer value iteration. The sample models Planner.hb, Solver.hb and Dyna.hb use this module.
Notation. Throughout, \(P \in [0, 1]^{S \times A \times S}\) (the argument p) has entries \(P_{sat}\) with \(\sum_t P_{sat} = 1\), \(r \in \mathbb{R}^{S \times A}\) has entries \(r_{sa}\), and \(\gamma \in [0, 1)\) (the argument g) is the discount. A policy \(\pi \in [0, 1]^{S \times A}\) (the argument pol) is row-stochastic.
open Prelude.Decisions
Shapes
Building blocks
stateOnes(p)— Let \(P \in \mathbb{R}^{S \times A \times S}\). Then \(\operatorname{stateOnes}(P) = \mathbf{1}_S\). A vector of ones, one per state, taken from the shape of p.flatOf(p)— Let \(P \in \mathbb{R}^{S \times A \times S}\). Then \(\operatorname{flatOf}(P)_{sa} = \sum_{t=1}^{S} 1 = S\) for every \(s, a\). A constant matrix of shape S×A, taken from the shape of p; every entry is \(S\) (not one), which is harmless where it is used, becauseuniformOfnormalises each row.
The backup and value iteration
The Bellman optimality backup is
\[Q(s,a) = r(s,a) + \gamma \sum_t p(t \mid s,a)\,V(t), \qquad V'(s) = \max_a Q(s,a),\]and value iteration repeats it from \(V = 0\).
backup
backup(p, r, g, v) · p transitions tensor S×A×S · r rewards matrix S×A · g discount scalar · v state values vector S → action values matrix S×A
Let \(P \in [0, 1]^{S \times A \times S}\) (the argument p), \(r \in \mathbb{R}^{S \times A}\), \(\gamma \in [0, 1)\) (the argument g) and \(v \in \mathbb{R}^{S}\) the state values. Then
One-step lookahead: the action values implied by the state values v.
Intuition “If I take this action now and the world is worth v afterwards, what is it worth?” Pay the immediate reward, then average the value of wherever you might land, discounted by one step. Every planning method on this page is this lookahead repeated.
Use it when you turn state values into action values: to act greedily on a value estimate, to plan inside Dyna with a learned model (modelled, payoffs), or as the building block of your own iteration. It costs \(S^2 A\) per call. For the value of the best action directly use sweep; for \(Q - V\) use advantage.
planningvalue-targets
def backup(p, r, g, v) ≔ r + g · einsum("sat,t->sa", p, v)
sweep
sweep(p, r, g, v) · p transitions tensor S×A×S · r rewards matrix S×A · g discount scalar · v state values vector S → state values vector S
Let \(P\), \(r\), \(\gamma\) and \(v \in \mathbb{R}^{S}\) be as in backup. Then
One sweep of value iteration: the maximum over actions of backup.
Intuition Each state asks every action “what are you worth if the future is v?” and keeps the best answer. Applied once, it makes the values correct for looking one step further ahead; the error to \(V^*\) shrinks by at least a factor \(\gamma\) per sweep.
Use it when you want control over the loop yourself: warm-starting from previous values after the model changed (a few Dyna sweeps per step), interleaving sweeps with learning, or checking convergence with drift. To run \(k\) sweeps from zero in one call use values.
planningcontrol
def sweep(p, r, g, v) ≔ rowMost(backup(p, r, g, v))
values
values(k, p, r, g) · k number of sweeps scalar · p transitions tensor S×A×S · r rewards matrix S×A · g discount scalar → optimal state values vector S
Let \(k \in \mathbb{N}\) be the number of sweeps and \(P\), \(r\), \(\gamma\) as above. Let \(V_0 = 0 \in \mathbb{R}^S\) and \(V_{j+1} = \operatorname{sweep}(P, r, \gamma, V_j)\) for \(j = 0, 1, \dots, k-1\). Then
\[\operatorname{values}(k, P, r, \gamma) = V_k \approx V^* .\]\(V^*\) approximated by \(k\) sweeps of value iteration from zero.
Intuition After \(k\) sweeps from zero, \(V_k(s)\) is the best total discounted reward achievable in the next \(k\) steps. As \(k\) grows the horizon passes the point where \(\gamma^k\) makes further rewards negligible, and \(V_k\) settles on \(V^*\).
Use it when you know (or have learned) \(P\) and \(r\) and want what each state is worth under optimal play. Pick \(k\) from the formula in the introduction — roughly horizon(g) from Returns times \(\ln(2 r_{\max} / (\varepsilon(1-\gamma)))\) — and confirm with drift. Cost is \(k \cdot S^2 A\). For very large \(\gamma\) it needs many sweeps but stays reliable; policyIteration needs fewer rounds, each a linear solve.
planningcontrolpolicy-evaluation
def values(k, p, r, g) ≔ valueLoop(k, p, r, g, zeroed(stateOnes(p)))
qValues
qValues(k, p, r, g) · k number of sweeps scalar · p transitions tensor S×A×S · r rewards matrix S×A · g discount scalar → optimal action values matrix S×A
Let \(k \in \mathbb{N}\), \(P\), \(r\), \(\gamma\) be as above and \(V_k = \operatorname{values}(k, P, r, \gamma)\). Then
\[\operatorname{qValues}(k, P, r, \gamma)_{sa} = r_{sa} + \gamma \sum_{t=1}^{S} P_{sat}\, (V_k)_t \approx Q^*_{sa} .\]\(Q^*\): one backup of the values from values.
Intuition State values say how good a place is; action values say how good each move from it is. One extra lookahead from \(V_k\) converts the first into the second, so the table ranks actions directly.
Use it when you need to compare actions, not just states: seeding a Q-table before learning, producing ground truth to test a Q-learner against, or softening the plan with boltzmann(qValues(k, p, r, g), tau). For just the greedy policy call plan; for the gap between each action and its state use advantage.
planningaction-selectionvalue-targets
def qValues(k, p, r, g) ≔ backup(p, r, g, values(k, p, r, g))
plan
plan(k, p, r, g) · k number of sweeps scalar · p transitions tensor S×A×S · r rewards matrix S×A · g discount scalar → policy matrix S×A
Let \(k \in \mathbb{N}\), \(P\), \(r\), \(\gamma\) be as above and \(Q = \operatorname{qValues}(k, P, r, \gamma)\). Then
\[\operatorname{plan}(k, P, r, \gamma)_{sa} = \operatorname{greedyOf}(Q)_{sa} = \frac{[\,Q_{sa} = \max_b Q_{sb}\,]}{\sum_{a'} [\,Q_{sa'} = \max_b Q_{sb}\,]} .\]The greedy policy with respect to qValues; tied actions share probability.
Intuition Once values are right, acting well is easy: in each state do whatever the one-step lookahead rates best. The policy is often correct long before the values have converged, because only the ranking of actions matters.
Use it when you want the answer to “what should I do in each state?” for a known or learned model — a planning agent, a Dyna step, or the reference policy for evaluating a learner. It returns a policy table, ready for evaluated or for rolling out with rollout. Near-ties can flip with too few sweeps, so check drift if the plan looks unstable.
planningcontroldecision-making
def plan(k, p, r, g) ≔ greedyOf(qValues(k, p, r, g))
drift
drift(p, r, g, v) · p transitions tensor S×A×S · r rewards matrix S×A · g discount scalar · v state values vector S → residual scalar
Let \(P\), \(r\), \(\gamma\) be as above and \(v \in \mathbb{R}^S\). Then
\[\operatorname{drift}(P, r, \gamma, v) = \max_{s} \big\lvert \operatorname{sweep}(P, r, \gamma, v)_s - v_s \big\rvert .\]The largest change one more sweep would make. Use it to check that \(k\) sweeps were enough: the values after one more sweep lie within \(\gamma/(1-\gamma)\) times this residual of \(V^*\).
Intuition A value table that is already optimal does not change under a sweep, so how much it would still change measures how far you are from done. Because each sweep contracts by \(\gamma\), a small residual guarantees a small distance to \(V^*\).
Use it when you chose a fixed sweep count and need to know whether it was enough, or you monitor a learned table’s consistency with a model. Compare it with your tolerance times \((1-\gamma)/\gamma\). It only checks the optimality equation, so it says nothing about a model that is itself wrong.
convergence-diagnosticsplanning
def drift(p, r, g, v) ≔ maximum(abs(sweep(p, r, g, v) − v))
Building blocks
valueLoop(k, p, r, g, v)— Let \(k \in \mathbb{N}\), \(P\), \(r\), \(\gamma\) as above and \(v \in \mathbb{R}^S\). Then \(\operatorname{valueLoop}(0, P, r, \gamma, v) = v\) and \(\operatorname{valueLoop}(k, P, r, \gamma, v) = \operatorname{valueLoop}(k - 1, P, r, \gamma, \operatorname{sweep}(P, r, \gamma, v))\) for \(k \ge 1\). k sweeps from the starting values v; used byvalues.
Exact policy evaluation
For a fixed policy \(\pi\) (S×A), the value solves the linear system
\[V^\pi = (I - \gamma P^\pi)^{-1}\, r^\pi, \qquad P^\pi_{st} = \sum_a \pi(a \mid s)\,p(t \mid s,a), \quad r^\pi_s = \sum_a \pi(a \mid s)\,r(s,a).\]underPolicy
underPolicy(p, pol) · p transitions tensor S×A×S · pol policy matrix S×A → state transition matrix matrix S×S
Let \(P \in [0, 1]^{S \times A \times S}\) (the argument p) and \(\pi \in [0, 1]^{S \times A}\) (the argument pol). Then
The Markov chain the policy induces on states.
Intuition Once the policy is fixed, the decisions disappear: from each state you move to \(t\) with the policy-weighted average of each action’s transition probabilities. What remains is a plain Markov chain over states.
Use it when you evaluate a fixed policy (it is the \(P^\pi\) in \(V^\pi\)), predict where the agent will be after \(n\) steps (multiply a state distribution by it \(n\) times), or find the long-run share of time spent in each state. Pair it with payUnder for the reward side.
policy-evaluationmarkov-chain-analysis
def underPolicy(p, pol) ≔ einsum("sa,sat->st", pol, p)
payUnder
payUnder(r, pol) · r rewards matrix S×A · pol policy matrix S×A → expected reward per state vector S
Let \(r \in \mathbb{R}^{S \times A}\) and \(\pi \in [0, 1]^{S \times A}\) (the argument pol). Then
The expected one-step reward in each state.
Intuition The reward half of collapsing a policy: in each state, what the agent collects on average next step, each action’s reward weighted by how often the policy picks it.
Use it when you evaluate a fixed policy (it is the \(r^\pi\) in \(V^\pi\)), compare the immediate payoff of two policies state by state, or compute the average reward rate against the chain’s long-run state distribution from underPolicy. It says nothing about future rewards; for those use evaluated.
policy-evaluationexpected-reward
def payUnder(r, pol) ≔ einsum("sa,sa->s", pol, r)
evaluated
evaluated(p, r, g, pol) · p transitions tensor S×A×S · r rewards matrix S×A · g discount scalar · pol policy matrix S×A → policy values vector S
Let \(P\), \(r\) and \(\gamma \in [0, 1)\) be as above and \(\pi \in [0, 1]^{S \times A}\) (the argument pol), and let \(M = I_S - \gamma P^{\pi}\). Then
\(V^\pi\) by one conjugate-gradient solve through the normal equations (see the note on conditioning above).
Intuition A policy’s value satisfies “my value is my expected reward plus \(\gamma\) times the value of where I go”. Written for all states at once that is a linear system, so instead of iterating toward \(V^\pi\) you solve for it in one shot.
Use it when you need the exact worth of a specific policy — scoring a hand-written rule, checking a learned policy against ground truth, or the evaluation half of improved. It costs one solve of an S×S system, fine for hundreds of states. Its accuracy degrades as \(\gamma\) nears 1 because the normal equations square the conditioning; past \(\gamma = 0.999\) iterate instead. To value the optimal policy use values.
policy-evaluationplanning
def evaluated(p, r, g, pol) ≔ solve(einsum("su,sv->uv", leftOver(p, g, pol), leftOver(p, g, pol)), einsum("su,s->u", leftOver(p, g, pol), payUnder(r, pol)))
Building blocks
leftOver(p, g, pol)— Let \(P \in [0, 1]^{S \times A \times S}\), \(\gamma \in [0, 1)\) and \(\pi \in [0, 1]^{S \times A}\). Then \(\operatorname{leftOver}(P, \gamma, \pi) = I_S - \gamma\, P^{\pi}\). The system matrix of policy evaluation.
Policy iteration
Policy iteration alternates exact evaluation with greedy improvement. Each round is a linear solve rather than a sweep, so a handful of rounds usually reaches the optimal policy.
improved
improved(p, r, g, pol) · p transitions tensor S×A×S · r rewards matrix S×A · g discount scalar · pol policy matrix S×A → policy matrix S×A
Let \(P\), \(r\), \(\gamma\) be as above, \(\pi \in [0, 1]^{S \times A}\) (the argument pol), and \(V^{\pi} = \operatorname{evaluated}(P, r, \gamma, \pi)\). Then
One improvement step: the greedy policy (ties split equally) with respect to the backup of \(V^\pi\).
Intuition Knowing exactly what your current policy is worth, ask in each state “is there a move that beats what I normally do, assuming I carry on as usual afterwards?” Switching to it can only help (the policy improvement theorem), and if nothing changes the policy is already optimal.
Use it when you improve a given policy one step at a time — for example to see how far a hand-designed rule is from optimal, or to run your own policy iteration with a custom stop test (stop when improved returns the same policy). Each call does a full evaluated solve. For \(k\) rounds from scratch use policyIteration.
controlplanningdecision-making
def improved(p, r, g, pol) ≔ greedyOf(backup(p, r, g, evaluated(p, r, g, pol)))
policyIteration
policyIteration(k, p, r, g) · k number of rounds scalar · p transitions tensor S×A×S · r rewards matrix S×A · g discount scalar → policy matrix S×A
Let \(k \in \mathbb{N}\) be the number of rounds and \(P\), \(r\), \(\gamma\) as above. Let \(\pi_0 = \tfrac{1}{A}\,\mathbf{1}_S \mathbf{1}_A^{\top}\) and \(\pi_{j+1} = \operatorname{improved}(P, r, \gamma, \pi_j)\) for \(j = 0, 1, \dots, k-1\). Then
\[\operatorname{policyIteration}(k, P, r, \gamma) = \pi_k .\]\(k\) rounds of policy iteration starting from the uniform random policy.
Intuition Evaluate, improve, repeat. A finite MDP has finitely many deterministic policies and every round is at least as good as the last, so the sequence reaches the optimum after a few rounds — usually far fewer rounds than value iteration needs sweeps.
Use it when you plan in a known finite MDP and prefer a handful of exact solves to hundreds of cheap sweeps; a single-digit \(k\) often suffices, and extra rounds only cost time. Each round costs a linear solve, so for large \(S\), or \(\gamma\) very close to 1 (see the conditioning note), plan with value iteration is the safer choice.
planningcontrol
def policyIteration(k, p, r, g) ≔ policyLoop(k, p, r, g, uniformOf(flatOf(p)))
Building blocks
policyLoop(k, p, r, g, pol)— Let \(k \in \mathbb{N}\), \(P\), \(r\), \(\gamma\) as above and \(\pi \in [0, 1]^{S \times A}\). Then \(\operatorname{policyLoop}(0, P, r, \gamma, \pi) = \pi\) and \(\operatorname{policyLoop}(k, P, r, \gamma, \pi) = \operatorname{policyLoop}(k - 1, P, r, \gamma, \operatorname{improved}(P, r, \gamma, \pi))\) for \(k \ge 1\). k improvement rounds from the policy pol; used bypolicyIteration.
Advantage
advantage
advantage(p, r, g, v) · p transitions tensor S×A×S · r rewards matrix S×A · g discount scalar · v state values vector S → advantages matrix S×A
Let \(P\), \(r\), \(\gamma\) be as above and \(v \in \mathbb{R}^S\) the state values. Then
\[\operatorname{advantage}(P, r, \gamma, v)_{sa} = r_{sa} + \gamma \sum_{t=1}^{S} P_{sat}\, v_t - v_s .\]What each action is worth over and above its state, \(A(s,a) = Q(s,a) - V(s)\), with \(Q\) the backup of v. At the optimal values the best action in each state has advantage 0 and the others are negative.
Intuition Subtracting the state’s own value removes “how lucky am I to be here” and leaves “how much better or worse is this move than the baseline”. With optimal values it reads as the cost of each mistake: \(-A(s, a)\) is what one wrong step loses before resuming optimal play.
Use it when you rank actions independently of state quality, measure the regret of a suboptimal action, find states where the choice barely matters (all advantages near 0), or compute exact advantages to test an estimator like gae against. It needs the model; from sampled trajectories use gae from Returns.
advantage-estimationdecision-makingregret-analysis
def advantage(p, r, g, v) ≔ backup(p, r, g, v) − lifted(backup(p, r, g, v), v)
Examples
The examples below live in one file, Planning.hb, so the class is Planning. They share one small MDP: two states and two actions, where action 0 stays and action 1 moves to the other state, and staying in state 1 pays 1. The best plan is to move from state 0 and then stay in state 1 forever, so \(V^*(1) = 1/(1-0.9) = 10\) and \(V^*(0) = 0.9 \cdot 10 = 9\).
using MLambda.Hilbert.Runtime;
var p = new double[2, 2, 2];
p[0, 0, 0] = 1; p[0, 1, 1] = 1; // from state 0: stay → 0, move → 1
p[1, 0, 1] = 1; p[1, 1, 0] = 1; // from state 1: stay → 1, move → 0
var pt = (Tensor)p;
var r = (Tensor)new double[,] { { 0, 0 }, { 1, 0 } };
Value iteration and the plan it implies
open Prelude.Decisions
fn optimalValues(k, p, r, g) ↦ values(k, p, r, g)
fn optimalPlan(k, p, r, g) ↦ plan(k, p, r, g)
fn sweepChange(p, r, g, v) ↦ drift(p, r, g, v)
using Examples;
var v = Planning.OptimalValues(200.0, pt, r, 0.9);
var plan = Planning.OptimalPlan(200.0, pt, r, 0.9);
var check = Planning.SweepChange(pt, r, 0.9, (Tensor)new[] { 9d, 10d });
Output (verified; the plan is shown flat in row order — state 0 moves, state 1 stays):
values = [9, 10]
plan = shape [2, 2] values [0, 1, 1, 0]
drift = 0
Evaluating the random policy exactly
Under the 50/50 policy the expected reward is \([0, 0.5]\) and both states lead anywhere with equal probability, so the average value is \(0.25/(1-0.9) = 2.5\) and the states differ by the immediate 0.5: \(V^\pi = [2.25, 2.75]\).
open Prelude.Decisions
fn exactValues(p, r, g, pol) ↦ evaluated(p, r, g, pol)
using Examples;
var random = (Tensor)new double[,] { { 0.5, 0.5 }, { 0.5, 0.5 } };
var vpi = Planning.ExactValues(pt, r, 0.9, random);
Output (verified):
evaluated = [2.25, 2.75]
Policy iteration and advantages
Three rounds of policy iteration reach the same plan as value iteration. At the optimal values, the chosen actions have advantage 0; staying in state 0 costs 0.9 and leaving state 1 costs 1.9.
open Prelude.Decisions
fn iteratedPolicy(k, p, r, g) ↦ policyIteration(k, p, r, g)
fn actionAdvantage(p, r, g, v) ↦ advantage(p, r, g, v)
using Examples;
var policy = Planning.IteratedPolicy(3.0, pt, r, 0.9);
var adv = Planning.ActionAdvantage(pt, r, 0.9, (Tensor)new[] { 9d, 10d });
Output (verified):
policyIteration = shape [2, 2] values [0, 1, 1, 0]
advantage = shape [2, 2] values [-0.9, 0, 0, -1.9]