Markov Chains & MDPs
The lesson has two parts. Part 1 is about chains: a transition matrix, the long-run distribution a chain settles into, PageRank, and the Metropolis sampler that turns a chain into a tool for drawing random samples. Part 2 adds decisions (Markov decision processes). Part 2 is optional if you are not heading toward reinforcement learning, and the Reinforcement Learning track picks it up from where this lesson stops.
After this lesson, you will be able to:
- Understand the Markov property — the future depends only on the present, not the path you took to get there — and recognise where it shows up in modern AI
- Compute transition matrices, k-step transitions, and stationary distributions, and see why solving for the stationary distribution is just an eigenvector problem from a previous lesson
- Compute PageRank on a tiny web (including a page with no outgoing links) and explain why the matrix is flipped to column form in the code
- Run a Metropolis sampler by hand on four states and in code on a two-bump target, and state the acceptance rule min(1, p(x') / p(x))
- Optional Part 2: solve a two-state decision process by hand with V = (I - gamma P)^-1 R, then read the general Bellman equations
- Connect Markov chains to MCMC, diffusion-model forward processes, and autoregressive LLM decoding without over-claiming what the Markov property alone gives you
Before You Start
#Part 1: Markov Chains
Part 1 needs only probability and one idea from the eigenvalues lesson. It ends with two payoffs: PageRank, and a sampler that turns a chain into a tool for drawing random numbers from hard distributions.
#The Present Screens Off the Past
This is the equation that runs reinforcement learning. The Markov property looks deceptively simple, but stacking it on top of itself produces some of the deepest equations in modern AI: the Bellman equation, the diffusion forward process, the autoregressive decoding loop in every LLM. Get this lesson and you can read RL, diffusion, and decoding papers without flinching.
X_0, X_1, X_2, … where each X_{t+1} depends only on X_t, not on the full history. Formally:Try it! Open the Python REPL (bottom-right of the screen: click Quick Actions, then Python) and type these lines yourself.
A chess engine is deciding its next move from the current board position. Does the engine need to know the full move history to play optimally, or is the position alone enough?
#A Worked Example: The Weather Chain
Suppose tomorrow's weather depends only on today's, with these transition probabilities:
- If today is sunny, tomorrow is sunny with probability 0.8 and rainy with probability 0.2
- If today is rainy, tomorrow is sunny with probability 0.4 and rainy with probability 0.6
P, where P_ij = "probability of going from state i to state j":π_0 = [1, 0] (definitely sunny), tomorrow's distribution is:Note on convention. This lesson uses row vectors (π_{t+1} = π_t P) because that is what most RL textbooks use. Some linear-algebra textbooks (and PyTorch code that prefers column vectors) write the same thing asπ_{t+1}^T = P^T π_t^T. They are mathematically identical — just different sides of the same product.
P, watch the state distribution π_t evolve one step at a time, and see how the particles redistribute. Try the "Sticky", "Mixing" and "Absorbing" presets to feel how the dynamics change.#Two Steps, Three Steps, k Steps
- Sunny → Sunny → Sunny: 0.8 × 0.8 = 0.64
- Sunny → Rainy → Sunny: 0.2 × 0.4 = 0.08
- Total: 0.72
P^2. The two-step transition matrix is P raised to the second power.P^2 = [[0.72, 0.28], [0.56, 0.44]]. So sunny → sunny in two steps has probability 0.72, matching our hand calculation. As k grows, something interesting happens to P^k: the rows start to look identical. Try it for k = 50:P50 = np.linalg.matrix_power(P, 50)
# [[0.667 0.333]
# [0.667 0.333]][0.667, 0.333]. The starting state has been forgotten. That repeated row is the stationary distribution.If P is a valid 3x3 transition matrix and you compute P^100, what structural property would you expect to see (assuming the chain is ergodic)?
#Stationary Distribution: The Long-Run Truth
π is a state distribution that is invariant under one step of the chain:π = π P says that π is a left eigenvector of P with eigenvalue 1. Solving for the stationary distribution is the same eigenvector problem you saw in the Eigenvalues & SVD lesson — there is no new linear algebra here, only a new application.P sums to 1, so multiplying P by a column of ones just returns that column: for the weather chain, P · [1, 1]ᵀ = [0.8 + 0.2, 0.4 + 0.6]ᵀ = [1, 1]ᵀ. So 1 is an eigenvalue of P, and Pᵀ has the same eigenvalues, so a left eigenvector for 1 exists. A result called the Perron-Frobenius theorem adds that a matrix with no negative entries has such an eigenvector with no negative entries either, and that is what lets us rescale it into probabilities.For the weather chain, working it out:
π = π Pgivesπ_S = 0.8 π_S + 0.4 π_Randπ_R = 0.2 π_S + 0.6 π_R- Both reduce to
π_S = 2 π_R - Combined with the normalisation
π_S + π_R = 1, we getπ_S = 2/3,π_R = 1/3
Over the long run, two-thirds of days are sunny, one-third rainy — regardless of where the chain started. That is the stationary distribution at work.
π_t oscillates instead of settling — the cyclic structure prevents convergence even though every state is visited equally often on average. Now try the "Random walk" preset and watch the same machinery quickly produce a uniform stationary distribution.P^t for several t and watches the rows align, then solves for the stationary distribution directly via eigendecomposition — the two answers must agree.#When Does π Exist? When Is It Unique?
π_t → π from any starting state:- Irreducibility: every state can reach every other state in finite time. There are no isolated islands.
- Aperiodicity: there is no fixed cycle length you are forced into. (Formally: the GCD of return times to any state is 1.)
π_t from any starting distribution. You do not need to memorise these conditions in detail — the takeaway is that "well-behaved" chains have a single long-run answer.π_t oscillates instead of converging, even though averages converge.P[s, s] is 1. A state is transient if the chain might leave it and never come back, and recurrent if the chain is certain to come back. In a chain with one absorbing state that every other state can reach, all the other states are transient, because the chain eventually falls into the absorbing state and stays.A 3-state chain has P[3, 3] = 1, so state 3 is absorbing. States 1 and 2 are not absorbing, and each of them can eventually reach state 3. What happens in the long run?
#Mixing Time (Briefly)
π_t converge to π? The answer is governed by the second-largest eigenvalue of P (in absolute value), often written |λ_2|. The gap 1 − |λ_2| is called the spectral gap, and convergence is exponentially fast at rate |λ_2|^t. A chain with a small spectral gap mixes slowly; a chain with a near-zero λ_2 snaps almost instantly to its stationary distribution. For the weather chain the eigenvalues of P are exactly 1 and 0.4, so |λ_2| = 0.4 and the gap is 0.6. The error should shrink like 0.4^t, and it does: after two steps the sunny-to-sunny entry is 0.72 against the limit 0.667, an error of 0.0533, which equals (1/3) × 0.4². After five steps the error is 0.0034.This single number controls how many MCMC samples you need, how quickly PageRank iterations converge, and how diffusion-model schedules trade off forward-process steps for sample quality.
#PageRank: The Web as a Markov Chain
Take four pages, A, B, C and D. A links to B and C. B links to C. C links to A and D. D links to nothing. The transition matrix, row by row (a row is "where do I go from this page"), starts like this:
| from \ to | A | B | C | D |
|---|---|---|---|---|
| A | 0 | 1/2 | 1/2 | 0 |
| B | 0 | 0 | 1 | 0 |
| C | 1/2 | 0 | 0 | 1/2 |
| D | ? | ? | ? | ? |
[1/4, 1/4, 1/4, 1/4].d = 0.85 the surfer follows a link, and with probability 1 − d = 0.15 the surfer teleports to a uniformly random page. Teleporting makes every page reachable from every other page, and it breaks any fixed cycle, so the chain is irreducible and aperiodic. That is exactly the ergodic condition from the last section, so the stationary distribution exists and is unique.i of P is where you go from state i, and the distribution updates as π P. PageRank code usually stores the transpose, M[j, i] = probability of going from page i to page j. Now each column of M sums to 1 (a column-stochastic matrix), and the update is π ← G π with the matrix on the left. It is the same chain written sideways: M is Pᵀ. Nothing changes except which way the table is read, but code that does it without saying so confuses everyone.M for these four pages, shows what goes wrong without the dangling-node repair, and then applies it.Run it and read it. Without the repair, the total probability collapses to about 2e-10, which means the "distribution" has drained away. With the repair, every column sums to 1 and the PageRank values are 0.234 for A, 0.187 for B, 0.345 for C and 0.234 for D, summing to 1, and the eigenvector gives the same four numbers. C ranks first because both A and B link to it. D ties with A: each of them receives half of C's clicks and a quarter of the jumps out of D. B ranks last: it receives only half of A's clicks and its quarter of the jumps out of D.
#The Metropolis Sampler: A Chain Built to Hit a Target
p(x) that you can evaluate but cannot sample from directly, for example a Bayesian posterior known only up to a constant. The trick is to build a Markov chain whose stationary distribution is p. Run the chain for a long time, and the fraction of time it spends near each x is p(x). The visited states are your samples. The Bayesian Inference in Practice lesson pointed ahead to this tool, and this is where it gets built.p as its stationary distribution? A sufficient condition is detailed balance. Write K(x → y) for the probability that the chain moves from x to y in one step. Detailed balance says that in equilibrium the flow each way is equal:K with this property using only ratios of p. From the current state x, propose a candidate x' from a symmetric proposal, meaning proposing x' from x is just as likely as proposing x from x'. Then accept the move with probabilityp(y) ≥ p(x). Going from x to y is always accepted, so K(x → y) = q, where q is the proposal probability. Going back from y to x is accepted with probability p(x)/p(y), so K(y → x) = q · p(x)/p(y). Then p(x) · q equals p(y) · q · p(x)/p(y), so the two flows match.1, 2, 4, 3. They add to 10, so p = [0.1, 0.2, 0.4, 0.3]. The proposal picks one of the other three states uniformly, which is symmetric. Start at state 2. Instead of random numbers, use these fixed proposals and fixed uniform draws u between 0 and 1. The rule is: accept when u is below the acceptance probability.| step | current x | proposed x' | ratio p(x')/p(x) | accept prob | u | decision | next state |
|---|---|---|---|---|---|---|---|
| 1 | 2 | 3 | 4/2 = 2 | 1 | 0.30 | accept | 3 |
| 2 | 3 | 1 | 1/4 = 0.25 | 0.25 | 0.80 | reject | 3 |
| 3 | 3 | 4 | 3/4 = 0.75 | 0.75 | 0.45 | accept | 4 |
| 4 | 4 | 2 | 2/3 = 0.667 | 0.667 | 0.10 | accept | 2 |
| 5 | 2 | 4 | 3/2 = 1.5 | 1 | 0.50 | accept | 4 |
| 6 | 4 | 1 | 1/3 = 0.333 | 0.333 | 0.60 | reject | 4 |
u is. Step 2 proposes a state four times less probable, so it is accepted only if u is below 0.25, and 0.80 is not. A rejection is not wasted: the chain stays where it is, and that stay counts as a visit. Over thousands of steps, state 3 (the heaviest) is visited about 40% of the time.p as its stationary distribution. Its transition matrix K has, for each row, K(x → y) = (1/3) · min(1, p(y)/p(x)) off the diagonal, with the leftover probability on the diagonal:K = [[0.0, 0.3333, 0.3333, 0.3333],
[0.1667, 0.1667, 0.3333, 0.3333],
[0.0833, 0.1667, 0.5, 0.25 ],
[0.1111, 0.2222, 0.3333, 0.3333]]
# p @ K = [0.1, 0.2, 0.4, 0.3] = p
# detailed balance for states 1 and 2:
# p(1) K(1 -> 2) = 0.1 * 1/3 = 0.0333
# p(2) K(2 -> 1) = 0.2 * 1/6 = 0.0333The same recipe works in one dimension with a continuous target. The cell below samples a two-bump distribution: 60% of the mass is a bump centred at -2 with spread 0.7, and 40% is a bump centred at +2 with spread 1.0. The proposal adds Gaussian noise, which is symmetric. The code compares the draws with the exact answer.
Read the output. About 48.4% of proposals were accepted. The sampled mean is -0.383 against an exact -0.400, and the sampled variance is 4.609 against an exact 4.534. The histogram sits within about 0.011 of the exact probability in every bin, including the narrow gap between the bumps: the chain does cross it, because a downhill step is sometimes accepted.
step to see how the acceptance rate and the effective sample size move together.#Hidden Markov Models, Briefly
X_t (the "true" state) and observations Y_t that are generated from X_t via an emission distribution. You see the Y_ts but not the X_ts; algorithms like Viterbi and Baum-Welch reconstruct the most likely hidden sequence or fit the parameters by EM. HMMs were the dominant speech-recognition technology before deep learning, and they remain useful in bioinformatics.#Part 2: Decisions and Markov Decision Processes (Optional)
Part 2 is optional if you are not heading toward reinforcement learning. Nothing in Part 1 depends on it, and the Reinforcement Learning track continues the story from here with Q-learning and the methods built on it. If you do read on, the plan is concrete first: after the definitions and the value functions, you solve a two-state decision problem by hand, and only then meet the general Bellman equations.
#Markov Decision Processes: Adding Actions and Rewards
a, the environment transitions to a new state, and the environment hands you a reward.(S, A, P, R, γ):- S: state space (where you can be)
- A: action space (what you can do)
- P(s' | s, a): transition function — the probability of landing in state
s'after taking actionain states - R(s, a) or R(s, a, s'): reward function — the immediate reward for the transition
- γ ∈ [0, 1]: discount factor — how much future reward is worth compared to immediate reward
π(a | s), which prescribes the probability of taking action a in state s. This π is not the same as the stationary-distribution π from earlier — yes, the symbol collides. RL textbooks lean on this overload heavily; just be aware.Why discount? Three reasons: (1) it keeps the sum finite even on infinite horizons; (2) it matches human / financial impatience; (3) it makes the math nice — the Bellman operator is a contraction, which we will see in a moment.
#Value Functions
#A Two-State MDP by Hand
P = [[0.8, 0.2], [0.4, 0.6]] with states Sunny and Rainy, and attach rewards: a sunny day is worth 1 point and a rainy day 0. To keep the arithmetic small, take γ = 0.5 and let the reward depend only on the state you are in, R = [1, 0]. With one action per state this is a Markov chain with rewards, and the value of each state is V(s), the expected discounted return from s.Write the value of Sunny in words: the point you get today, plus half of the average value of tomorrow. Tomorrow is Sunny with probability 0.8 and Rainy with 0.2, so:
V(S) = 1 + 0.5 × (0.8 V(S) + 0.2 V(R))V(R) = 0 + 0.5 × (0.4 V(S) + 0.6 V(R))
V = R + γ P V. Move the P V term across: (I − γ P) V = R, soI − γP = [[1 − 0.4, −0.1], [−0.2, 1 − 0.3]] = [[0.6, −0.1], [−0.2, 0.7]]. For a 2 by 2 matrix, the inverse swaps the diagonal entries, negates the off-diagonal ones, and divides by the determinant. The determinant is 0.6 × 0.7 − (−0.1)(−0.2) = 0.42 − 0.02 = 0.40. So the inverse is (1 / 0.4) × [[0.7, 0.1], [0.2, 0.6]] = [[1.75, 0.25], [0.5, 1.5]], and multiplying by R = [1, 0] picks out the first column:V = [1.75, 0.5]1 + 0.5 × (0.8 × 1.75 + 0.2 × 0.5) = 1 + 0.5 × 1.5 = 1.75, and 0 + 0.5 × (0.4 × 1.75 + 0.6 × 0.5) = 0.5 × 1.0 = 0.5. Both hold. The same thing in numpy:import numpy as np
P = np.array([[0.8, 0.2], [0.4, 0.6]])
R = np.array([1.0, 0.0])
gamma = 0.5
V = np.linalg.solve(np.eye(2) - gamma * P, R)
print(V) # [1.75 0.5 ]
print(R + gamma * P @ V) # [1.75 0.5 ] the equation holdsR + γ P V, is one Bellman backup. Start from V = [0, 0] and apply it repeatedly. The first backup gives [1, 0]. The second gives V(S) = 1 + 0.5 × (0.8 × 1 + 0.2 × 0) = 1.4 and V(R) = 0 + 0.5 × (0.4 × 1 + 0.6 × 0) = 0.2, so [1.4, 0.2]. Then [1.58, 0.34], [1.666, 0.418], [1.7082, 0.4586], creeping up on [1.75, 0.5]. Repeated backups and the direct solve agree.[1.75, 0.5] as a guide, compare the two actions in Rainy with one backup each:- Stay in:
0 + 0.5 × (0.4 × 1.75 + 0.6 × 0.5) = 0.5 - Go out:
−0.2 + 0.5 × (0.9 × 1.75 + 0.1 × 0.5) = −0.2 + 0.8125 = 0.6125
V* = [1.7714, 0.6286], and the best action in Rainy is still "go out". In Sunny there is only one action, so nothing to choose.#Bellman Equations: The Crown Jewel
s equals the expected immediate reward plus the discounted expected return from the next state:Q^π:V^* and a corresponding optimal policy π^* that achieves it. The reason it has a unique solution is that the Bellman optimality operator is a contraction.You are designing an MDP for a self-driving car. Episodes are not guaranteed to terminate (the road can be infinitely long). If you set the discount factor γ = 1.0, what is the most likely failure mode?
#Value Iteration and Policy Iteration
Two classical algorithms turn the Bellman equations into computational procedures:
- Value iteration: start with any
V_0. Repeatedly applyV_{k+1}(s) = max_a Σ_{s'} P(s'|s,a)[R + γ V_k(s')]. Converges toV^*geometrically. Extract the optimal policyπ^*(s) = argmax_a [...]at the end. - Policy iteration: alternate between (a) policy evaluation — solve the Bellman expectation equation to get
V^πfor the current π, and (b) policy improvement, updateπ(s) ← argmax_a Q^π(s, a). Converges in finitely many iterations on finite MDPs.
Q_θ. PPO is approximate policy iteration with a stochastic policy and a trust-region constraint. Knowing the Bellman equations means recognising every paper as a variation on this theme.Press Step to apply one Bellman sweep and watch the value spread outward from the goal one cell at a time, then raise the slip probability or lower gamma and watch the greedy arrows turn away from the trap.
#Connecting It Back to Modern AI
Put together everything in this lesson and you can read most of modern AI's "stochastic" math in one breath.
π(a|s) ∝ exp(Q(s,a) / τ), a softmax over the Q-values divided by a temperature τ (the Greek letter tau; a large τ makes the choice nearly uniform, a small one nearly greedy). A policy network in policy-gradient methods does something else. It outputs raw scores called logits, one per action, and F.softmax turns those logits into π(a | s). The logits are learned directly. They are not Q-values.q(x_t | x_{t-1}) = N(x_t; √(1−β_t) x_{t-1}, β_t I) is a Markov chain on images: each step shrinks the previous image a little and adds fresh Gaussian noise of size β_t (a small number chosen per step). The Markov property is what lets you describe the whole process as a composition of these one-step moves. The closed form for q(x_t | x_0) comes from something else, the linearity of Gaussians (multivariate-Gaussian lesson): a linear map of a Gaussian, plus independent Gaussian noise, is again Gaussian. Applying that 1000 times collapses the steps into one Gaussian, which is what makes the DDPM training objective a single-step regression and not a 1000-step Monte Carlo estimate.The bridges, stated explicitly so you can come back and find them:
- PageRank is the stationary distribution of the link Markov chain on the web. Google's original PageRank was an eigenvector computation (Eigenvalues & SVD lesson) on a stochastic transition matrix (this lesson), with the dangling-node and damping repairs from the PageRank section.
- Diffusion models' forward process
q(x_t | x_{t-1}) = N(x_t; √(1−β_t) x_{t-1}, β_t I)is a Markov chain of Gaussian steps. The chain structure lets you compose the steps, and the closed-formq(x_t | x_0)follows because linear Gaussian updates keep everything Gaussian. - Q-learning trains a neural network
Q_θ(s, a)to satisfy the Bellman optimality equation:Q_θ(s, a) ← R + γ max_{a'} Q_θ(s', a'). The Bellman optimality from this lesson is the fixed-point that DQN, Rainbow, and all value-based RL methods chase. - PPO and TRPO use trust-region constraints (Lagrangian / KL-bounded updates). The discount factor γ in the Bellman equation reappears in the GAE advantage estimator:
Â_t = Σ_l (γλ)^l δ_{t+l}, where δ is the TD error. Both terms are glossed in the DeepDive after this list. - MCMC samplers construct ergodic Markov chains whose stationary distribution is the target posterior. You built the simplest one, Metropolis with a symmetric proposal, in the Metropolis section; Metropolis-Hastings adds a correction for non-symmetric proposals, and Gibbs and Hamiltonian Monte Carlo are other ways to build such chains. Bayesian deep learning's posterior sampling, used in tools like PyMC and NumPyro, is built on these Markov-chain foundations.
Here is the PageRank from scratch again, on a bigger 5-page web with no dead ends. It is the same stationary-distribution computation as in the PageRank section, in about thirty lines, with the same column-stochastic convention.
Reinforcement Learning: An Introduction (2nd edition)
Richard S. Sutton, Andrew G. Barto (2018)
#Try It Yourself
You have met the stationary distribution three ways in prose: as the long-run share of time, as a left eigenvector, and as the limit of P^t. Here you check that they really agree on a fresh 3-state weather chain, then bolt on actions and rewards and run value iteration yourself. Work through the TODOs in order; the solution is there when you want to compare.
Tests · Verify the three stationary estimates agree to about 3 decimals: simulated [0.456, 0.282, 0.262], eigenvector and P^100 both [0.457, 0.283, 0.261]. For the stretch MDP, verify V = [16.317, 14.858, 15.23] and policy [0, 0, 1].
Read your output against these numbers. The simulation gives visit frequencies of 0.456, 0.282 and 0.262, while the eigenvector and the matrix power both give 0.457, 0.283 and 0.261. The simulation is off by at most 0.0014, which is ordinary Monte Carlo noise for 100,000 steps, and the other two methods agree exactly because they are the same fact computed differently. For the stretch, 50 sweeps with gamma 0.9 give V = [16.317, 14.858, 15.23], and the greedy policy is [0, 0, 1]: keep action 0 in Sunny and Cloudy, and switch to action 1 in Rainy. Try changing gamma to 0.5 and see whether the policy changes.
#Quick Check
What does the Markov property say about a process?