30 - Markov Chains: The Property That Forgets the Past, and the Distribution That Stops Changing

Machine Learning Engineer course · Module 02 · Mathematics and statistics · Lesson 15

Goal. Recognise when a process is a Markov chain, write down its transition matrix and compute its stationary distribution in three independent ways. By the end you can say whether a chain has a limiting distribution, compute it as a left eigenvector of eigenvalue 1, estimate how many steps it needs to settle by looking at the second eigenvalue, and tell the difference between a process you undergo and a process in which you decide.

A myth to dispel. «Markov means the past does not matter.» No: the Markov property says the past does not matter in addition to the present, because the state — if it is chosen well — already summarises all of it. It is not a property of the phenomenon: it is a property of the phenomenon-state pair. Choosing the state badly is the most common way of building a broken Markov model.

Introduction

This is the last lesson of the mathematics module, and it closes two circles at once. The first: in the lesson on eigenvalues, the SVD and PCA you learned that some matrices have special directions, and there the eigenvectors were used to compress data. Here they do something entirely different: the eigenvector of eigenvalue 1 of a transition matrix is the distribution a process tends to, and the second eigenvalue says how long it takes to get there. It is the only serious use of eigenvalues in this course outside PCA.

The second circle looks forward. In module 08 the Markov decision process is the language in which every reinforcement learning problem is written, and the Bellman equation is its consequence. Neither of those lessons, however, defines the chain underneath: the transition matrix, the stationary distribution, the condition for it to exist. We build it here, so that there only the decision has to be added.

Final result of this lesson. Given a graph of six pages with declared links, you can build the transition matrix, compute the stationary distribution in three ways that agree to 1e-10, estimate the settling time from the spectral gap, diagnose a page with no outgoing links and fix it with the teleportation factor.

Lesson roadmap.

1. The Markov property, and the processes that break it until you widen the state.
2. The transition matrix, the row-column convention and the six-page crawl graph.
3. The distribution after n steps, Chapman-Kolmogorov and computing Pⁿ by repeated squaring.
4. Communication classes, recurrent and transient states.
5. Irreducibility, period and aperiodicity, with the chain that oscillates forever.
6. The stationary distribution as a left eigenvector, Perron-Frobenius, three ways to compute it.
7. The spectral gap, convergence speed and mixing time, measured.
8. Detailed balance, and why chains are designed on top of it.
9. The dangling page problem and the teleportation factor.
10. From chain to decision process: the line that leads to module 08.

The outline has two halves: the first five entries build the vocabulary (property, matrix, powers, classes, period), the last five use it to answer the practical questions (what the stationary distribution is, how long it takes, how to design a chain, what happens if it breaks). If you are only interested in reinforcement learning, you can read the lesson up to the stationary distribution and then skip to the last entry.

The Markov property

The analogy: snakes and ladders

You are on square 34 of a snakes and ladders board. To know where you will land on the next throw you need one thing only: that you are on 34. It does not matter whether you got there with a 6 and a 2, with three unlucky throws or by a shortcut: the board does not remember, and neither does the die. The future depends on the past only through the current position.

A process of this kind is called a Markov chain, after the Russian mathematician Andrey Andreyevich Markov, who introduced them in 1906 to extend the law of large numbers to variables that depend on one another. It is the simplest structure that is not independence: not «every step is unrelated to the others», but «every step is related only to the previous one».

The formal definition

Definition. A sequence of random variables X₀, X₁, X₂, … taking values in a finite set of states is a Markov chain if, for every n and every history,

P(Xₙ₊₁ = j | Xₙ = i, Xₙ₋₁ = iₙ₋₁, …, X₀ = i₀) = P(Xₙ₊₁ = j | Xₙ = i).

Read it as a sentence: the probability of ending up in j at the next step, knowing the whole history, is the same as the one you get knowing only where you are now. If that probability does not depend on n either — that is, the rules of the game do not change over time — the chain is called homogeneous, and we write Pᵢⱼ = P(Xₙ₊₁ = j | Xₙ = i). Every chain in this lesson is homogeneous, and it is worth declaring it: it is an assumption, not a fact of nature.

The warning that prevents the classic mistake. The property does not say the past is irrelevant. It says the past is already summarised in the state. These are two very different claims: the first is almost always false in real systems, the second is a design requirement that it is up to you to satisfy, by choosing what goes inside the state.

Two processes that are not Markovian (until you widen the state)

First example: browsing with the «back» button. A visitor moves around a website and now and then presses «back», which returns them to the page they came from. If you take «the page I am on» as the state, the process is not Markovian: the probability of landing on /blog depends on where you were before, and that information is not in the state. The fix is not to abandon Markov, it is to widen the state: take the pair (current page, previous page) as the state. On pairs the property holds again, at the cost of going from 6 states to (at most) 36.

Second example: the sum of the last two die rolls. Let Sₙ = Dₙ + Dₙ₋₁. If Sₙ = 7, the distribution of Sₙ₊₁ = Dₙ₊₁ + Dₙ depends on the value of Dₙ: with Dₙ = 6 the next sum lies between 7 and 12, with Dₙ = 1 between 2 and 7. The state «sum» throws away precisely the piece of the past that is needed. Fix: use Dₙ as the state (or the pair (Dₙ₋₁, Dₙ) if you want to keep reading the sum off the state).

The rule to take away. «It is not Markovian» almost always means «the state is too poor». Widening the state has a price: the number of states explodes, and the transition matrix grows with the square. Choosing the state is the real design work of a Markov model, and it is paid for in memory.

The transition matrix

With a finite number of states, all one-step probabilities fit in a square matrix P, where Pᵢⱼ is the probability of going from i to j. Two properties define it: Pᵢⱼ ≥ 0 and, for every row,

Σⱼ Pᵢⱼ = 1.

In words: from any state, you end up somewhere. A matrix with these two properties is called row-stochastic. It is not a universal convention: half of the physics and linear algebra literature uses column-stochastic matrices, with Pᵢⱼ equal to the probability of going from j to i.

The convention, declared once and for all. In this lesson P is row-stochastic, and the probability distribution is a row vector: if π₀ is the initial distribution, after one step the distribution is π₁ = π₀ P, and after n steps

πₙ = π₀ Pⁿ.

Writing P @ pi instead of pi @ P is mistake number one for anyone implementing a chain for the first time: it raises no exception, because the matrix is square, and it simply returns the numbers of a different chain (the transposed one). Declare the convention in a comment above the matrix and check it immediately, verifying that the rows sum to 1 and not the columns.

The six-page crawl graph

The running example of this lesson is a mini-site of six pages, with the links declared by hand: no external dataset, no downloaded network. A crawler starts from a page and at every step follows an outgoing link chosen at random with uniform probability. The six pages are /, /blog, /blog/a1, /blog/a2, /strumenti, /contatti, and the link structure is the one in the figure.

Translating it into code is mechanical: a dictionary listing the outgoing links of each page, and a function that turns it into a matrix by spreading 1/k over the k links.

import numpy as np

pagine = ["/", "/blog", "/blog/a1", "/blog/a2", "/strumenti", "/contatti"]
link = {
    "/":           ["/blog", "/strumenti", "/contatti"],
    "/blog":       ["/blog/a1", "/blog/a2", "/"],
    "/blog/a1":    ["/blog", "/blog/a2"],
    "/blog/a2":    ["/blog", "/blog/a1"],
    "/strumenti":  ["/", "/contatti"],
    "/contatti":   ["/"],
}

def matrice_di_transizione(pagine, link):
    """P[i, j] = probability of moving from page i to page j in one click."""
    indice = {p: i for i, p in enumerate(pagine)}
    P = np.zeros((len(pagine), len(pagine)))
    for sorgente, uscite in link.items():
        for destinazione in uscite:
            P[indice[sorgente], indice[destinazione]] = 1.0 / len(uscite)
    return P

P = matrice_di_transizione(pagine, link)
for nome, riga in zip(pagine, P):
    print(f"{nome:>11} | " + "  ".join(f"{v:.3f}" for v in riga))
print("every row sums to 1:", bool(np.allclose(P.sum(axis=1), 1.0)))

Output obtained (Python 3.12, NumPy 2.3):

          / | 0.000  0.333  0.000  0.000  0.333  0.333
      /blog | 0.333  0.000  0.333  0.333  0.000  0.000
   /blog/a1 | 0.000  0.500  0.000  0.500  0.000  0.000
   /blog/a2 | 0.000  0.500  0.500  0.000  0.000  0.000
 /strumenti | 0.500  0.000  0.000  0.000  0.000  0.500
  /contatti | 1.000  0.000  0.000  0.000  0.000  0.000
every row sums to 1: True

What the code does. The link dictionary is the only thing written by hand: everything else is derived. matrice_di_transizione puts 1/len(uscite) in cell (source, destination), which is the assumption «the crawler picks a link at random, all with the same probability» — an assumption, not a fact: a real crawler follows links with different probabilities, and a human visitor even more so. The /contatti row has a single 1.000 because that page has exactly one outgoing link. The final check on the rows is the only defence against a typo in the dictionary.

The distribution after n steps

If today's distribution is the row vector π₀, tomorrow's is π₀P. Row i of P says where the mass sitting in i goes, and the row-by-matrix product adds up the contributions of all departures. Repeating, after n steps you get πₙ = π₀Pⁿ: the n-th power of the matrix holds all the n-step transition probabilities, because (Pⁿ)ᵢⱼ = P(Xₙ = j | X₀ = i).

def distribuzione_dopo(P, pi0, n):
    """Row vector times matrix, n times: π₀ Pⁿ."""
    pi = pi0.astype(float).copy()
    for _ in range(n):
        pi = pi @ P                      # row times matrix: NOT P @ pi
    return pi

pi0 = np.zeros(6)
pi0[pagine.index("/contatti")] = 1.0    # the crawler starts from /contatti

for n in (0, 1, 2, 3, 4, 100):
    pi = distribuzione_dopo(P, pi0, n)
    print(f"n = {n:>3} | " + "  ".join(f"{v:.4f}" for v in pi))

Output obtained:

n =   0 | 0.0000  0.0000  0.0000  0.0000  0.0000  1.0000
n =   1 | 1.0000  0.0000  0.0000  0.0000  0.0000  0.0000
n =   2 | 0.0000  0.3333  0.0000  0.0000  0.3333  0.3333
n =   3 | 0.6111  0.0000  0.1111  0.1111  0.0000  0.1667
n =   4 | 0.1667  0.3148  0.0556  0.0556  0.2037  0.2037
n = 100 | 0.2400  0.2400  0.1600  0.1600  0.0800  0.1200

What the code does. It starts from the certainty of being on /contatti and applies the matrix once per step. The first steps are deterministic or nearly so: from /contatti you must go to /, and from there you go with one third each to /blog, /strumenti, /contatti. After four steps the memory of the starting point has almost vanished, and after a hundred the distribution no longer moves by more than 1e-13: that is the stationary distribution. Careful: «it no longer moves» is not automatic, it is a consequence of two properties of the chain (irreducibility and aperiodicity) that we will define shortly and that, when missing, break everything.

Chapman-Kolmogorov, and computing Pⁿ

Going from i to j in m + n steps means passing through some intermediate state k after m steps. Summing over all possible k gives the Chapman-Kolmogorov equation:

(P^(m+n))ᵢⱼ = Σₖ (Pᵐ)ᵢₖ · (Pⁿ)ₖⱼ,

which is, word for word, the definition of matrix product: P^(m+n) = Pᵐ Pⁿ. It is not a deep theorem, it is the Markov property written in matrix notation. It does have a practical consequence though: computing P¹⁰⁰ does not take 99 multiplications, about ten are enough if you use exponentiation by repeated squaring, which splits the exponent into powers of two.

def potenza_veloce(M, n):
    """Mⁿ by repeated squaring; also returns how many multiplications it took."""
    R = np.eye(len(M))
    B = M.copy()
    moltiplicazioni = 0
    while n > 0:
        if n % 2 == 1:
            R = R @ B
            moltiplicazioni += 1
        n //= 2
        if n > 0:
            B = B @ B
            moltiplicazioni += 1
    return R, moltiplicazioni

P3, _ = potenza_veloce(P, 3)
P4, _ = potenza_veloce(P, 4)
P7, _ = potenza_veloce(P, 7)
print("Chapman-Kolmogorov, P⁷ = P³ P⁴ :", bool(np.allclose(P7, P3 @ P4)))
print("P⁷ = P⁴ P³ (same result)       :", bool(np.allclose(P7, P4 @ P3)))
P100, conta = potenza_veloce(P, 100)
print(f"P¹⁰⁰: {conta} multiplications instead of 99")
print("all the rows of P¹⁰⁰ are equal:", bool(np.allclose(P100, np.ones((6, 1)) @ P100[:1])))

Output obtained:

Chapman-Kolmogorov, P⁷ = P³ P⁴ : True
P⁷ = P⁴ P³ (same result)       : True
P¹⁰⁰: 9 multiplications instead of 99
all the rows of P¹⁰⁰ are equal: True

What it tells us. The last line is the most important of the four: in P¹⁰⁰ all rows coincide, that is, after a hundred steps it no longer matters where you started. A matrix of this kind has rank 1 in the limit, and the repeated row is exactly the stationary distribution. Trade-off. Repeated squaring pays off for a large n on a small dense matrix; if the matrix is huge and sparse (a real web graph) the powers fill it with zeros that turn into non-zeros, and it is better to multiply a single vector repeatedly, which costs one sparse product per step and never builds Pⁿ.

Classifying the states

Before asking «where does it end up», you need to know what the chain's graph looks like. Four definitions are enough. State j is accessible from i (written i → j) if there is an n ≥ 0 with (Pⁿ)ᵢⱼ > 0, that is, if the graph lets you go from i to j following the arrows. Two states communicate if i → j and j → i. Communication is an equivalence relation, so it splits the states into communication classes. A class is closed if no arrow leaves it.

From here follow the two labels that matter. A state is recurrent if, starting from it, the probability of coming back sooner or later is 1; it is transient if that probability is less than 1, and in that case the chain returns to it only finitely many times, then abandons it forever. On a finite state space there is a handy shortcut: the states of a closed class are recurrent, those of a non-closed class are transient. The chain is irreducible if there is a single class, that is, if from every page you can reach every other page.

def raggiungibili(P, i):
    """States reachable from i in any number of steps (graph traversal)."""
    visti, da_visitare = {i}, [i]
    while da_visitare:
        k = da_visitare.pop()
        for j in np.flatnonzero(P[k] > 0):
            if int(j) not in visti:
                visti.add(int(j))
                da_visitare.append(int(j))
    return visti

def classi_di_comunicazione(P):
    avanti = [raggiungibili(P, i) for i in range(len(P))]
    classi, assegnati = [], set()
    for i in range(len(P)):
        if i in assegnati:
            continue
        classe = sorted(j for j in range(len(P)) if j in avanti[i] and i in avanti[j])
        classi.append(classe)
        assegnati.update(classe)
    return classi

def chiusa(P, classe):
    """A class is closed if no transition leaves it: on finite states, it is recurrent."""
    return all(P[i, j] == 0 for i in classe for j in range(len(P)) if j not in classe)

def descrivi(etichetta, M):
    classi = classi_di_comunicazione(M)
    print(f"{etichetta}: {len(classi)} class(es), irreducible = {len(classi) == 1}")
    for classe in classi:
        tipo = "recurrent" if chiusa(M, classe) else "transient"
        print(f"   {tipo:>11} | {[pagine[i] for i in classe]}")

descrivi("crawl graph", P)

Output obtained:

crawl graph: 1 class(es), irreducible = True
     recurrent | ['/', '/blog', '/blog/a1', '/blog/a2', '/strumenti', '/contatti']

What the code does. raggiungibili is a graph traversal that ignores the values of the probabilities and only looks at which ones are non-zero: for accessibility the structure matters, not the numbers. classi_di_comunicazione groups the mutually reachable states; chiusa checks that no arrows leave. Our graph has a single class and is therefore irreducible: every page is recurrent, and the crawler will pass through each of them infinitely often. We will reuse this same descrivi function later on the broken version of the site, where the answer changes.

Irreducibility, period and aperiodicity

Irreducibility alone is not enough to guarantee that the distribution settles. A second property is needed. The period of a state i is

d(i) = gcd{ n ≥ 1 : (Pⁿ)ᵢᵢ > 0 },

that is, the greatest common divisor of the lengths of the paths returning to i. If d(i) = 1 the state is called aperiodic. In an irreducible chain all states have the same period, so one can speak of the period of the chain. Any self-loop (Pᵢᵢ > 0) makes the chain aperiodic at once, because it puts the number 1 into the set of return times; but it is not necessary: two cycles of coprime lengths are enough, like our / → /strumenti → / of length 2 and / → /strumenti → /contatti → / of length 3.

from math import gcd

def periodo(P, i, massimo=24):
    """gcd of the lengths n for which it is possible to return to i in n steps."""
    d, Pn = 0, np.eye(len(P))
    for n in range(1, massimo + 1):
        Pn = Pn @ P
        if Pn[i, i] > 1e-12:
            d = gcd(d, n)
    return d

print("period of each page of the crawl graph:", [periodo(P, i) for i in range(6)])

# two pages linking only to each other: a periodic chain of period 2
Q = np.array([[0.0, 1.0],
              [1.0, 0.0]])
print("period of the two pages:", [periodo(Q, i) for i in range(2)])
pi = np.array([1.0, 0.0])
for n in range(5):
    print(f"   n = {n}: {pi[0]:.2f} {pi[1]:.2f}")
    pi = pi @ Q
pi_stazionaria = np.array([0.5, 0.5])
print("π = (0.5, 0.5) satisfies π = πQ:", bool(np.allclose(pi_stazionaria @ Q, pi_stazionaria)))
print("moduli of the eigenvalues:", "  ".join(f"{m:.4f}" for m in sorted(abs(np.linalg.eigvals(Q)), reverse=True)))

Output obtained:

period of each page of the crawl graph: [1, 1, 1, 1, 1, 1]
period of the two pages: [2, 2]
   n = 0: 1.00 0.00
   n = 1: 0.00 1.00
   n = 2: 1.00 0.00
   n = 3: 0.00 1.00
   n = 4: 1.00 0.00
π = (0.5, 0.5) satisfies π = πQ: True
moduli of the eigenvalues: 1.0000  1.0000

What it tells us. The two-page chain is irreducible — you go from one to the other and back — and yet the distribution never converges: it oscillates between (1, 0) and (0, 1) forever. The vector (0.5, 0.5) is stationary, in the sense that it satisfies π = πQ, but no initial distribution other than itself ever reaches it. The numerical tell is the last line: there are two eigenvalues of modulus 1 (they are 1 and −1), and that −1 eigenvalue is exactly the oscillation. The rule to take away. Irreducible guarantees that the stationary distribution exists and is unique; aperiodic guarantees that you get there. They are two separate requirements, and the second is the one more often forgotten.

The stationary distribution

Definition. A probability distribution π (row vector, non-negative components summing to 1) is stationary if

π = π P.

In words: if today the probability of being on each page is π, tomorrow it will still be π. It does not mean the crawler stands still — it keeps jumping from one page to another — it means that the inflow of each page exactly balances its outflow. It is a dynamic equilibrium, like the level of a bath with both the tap and the drain open.

It is a left eigenvector, and this closes the circle with eigenvalues

Transposing both sides of π = πP gives

Pᵀ πᵀ = πᵀ,

that is: πᵀ is an eigenvector of Pᵀ for the eigenvalue 1. Equivalently, π is a left eigenvector of P with eigenvalue 1. Eigenvalue 1 is always there, for a one-line reason: the column vector of all ones 𝟙 satisfies P𝟙 = 𝟙 (that is the rows-sum-to-1 condition), so 1 is an eigenvalue of P, and a matrix and its transpose have the same eigenvalues. The vector 𝟙 is the right eigenvector; the left one, which looks nothing like it, is the stationary distribution.

Here you can see why in the lesson on eigenvalues we insisted on the distinction between right and left eigenvectors. In PCA the matrix was symmetric and the distinction disappeared: for a symmetric matrix right and left coincide. A transition matrix, instead, is almost never symmetric — /contatti points to / with probability 1, / points to /contatti with probability one third — and swapping the two eigenvectors is a mistake that produces plausible, wrong numbers.

Perron-Frobenius, stated precisely

Theorem (Perron 1907, Frobenius 1912), stochastic matrix case. Let P be the transition matrix of a finite, irreducible Markov chain. Then: (a) there exists one and only one probability distribution π with π = πP, and all its components are strictly positive; (b) every eigenvalue has modulus ≤ 1, and the eigenvalue 1 is simple, that is, of multiplicity 1; (c) if the chain is also aperiodic, then 1 is the only eigenvalue of modulus 1, and πₙ → π for any initial distribution π₀.

We state it without proving it: the standard proof is in Norris, Markov Chains, chapter 1, and the general version for non-negative matrices in Horn and Johnson, Matrix Analysis, chapter 8. We will check it numerically, which is a different thing and is worth less than a proof, but it serves to confirm that you have understood the statement. It is worth noting where each hypothesis enters: without irreducibility the stationary distribution may not be unique (two separate closed classes give infinitely many, all their combinations); without aperiodicity it exists and is unique, but you never reach it, as in the two-page example.

Three ways to compute it

The first is the most direct: π = πP is a homogeneous linear system π(P − I) = 0, which on its own has infinitely many solutions (all multiples of π). It is made unique by replacing any one of the equations — which is redundant — with the constraint Σπᵢ = 1. The second asks a library for the eigenvector. The third multiplies a vector by the matrix until it stops moving.

def stazionaria_sistema_lineare(P):
    """Solves (Pᵀ − I)πᵀ = 0, replacing the last equation with Σπᵢ = 1."""
    n = len(P)
    A = P.T - np.eye(n)
    A[-1, :] = 1.0                      # the redundant equation becomes the normalisation
    b = np.zeros(n)
    b[-1] = 1.0
    return np.linalg.solve(A, b)

def stazionaria_autovettore(P, mostra=False):
    """Left eigenvector of P for eigenvalue 1 = right eigenvector of Pᵀ."""
    valori, vettori = np.linalg.eig(P.T)                 # careful: COMPLEX numbers
    if mostra:
        print("eigenvalues of Pᵀ (complex):",
              "  ".join(f"{v.real:+.4f}{v.imag:+.4f}i" for v in valori))
    k = int(np.argmin(np.abs(valori - 1.0)))             # the closest to 1, not "the first"
    if mostra:
        print(f"chosen index: {k} | eigenvalue: {valori[k].real:+.6f}{valori[k].imag:+.6f}i")
        print("raw eigenvector (real part):", "  ".join(f"{v:+.4f}" for v in np.real(vettori[:, k])))
    v = np.real(vettori[:, k])                           # drop the imaginary part, which is ~0
    return v / v.sum()                  # without this line the norm is 1, not the sum

def stazionaria_potenze(P, passi=1000):
    """Power method on the row vector, starting from the uniform distribution."""
    pi = np.full(len(P), 1.0 / len(P))
    for _ in range(passi):
        pi = pi @ P
    return pi

metodi = {
    "linear system": stazionaria_sistema_lineare(P),
    "left eigenvector": stazionaria_autovettore(P, mostra=True),
    "power method": stazionaria_potenze(P),
}
for nome, pi in metodi.items():
    print(f"{nome:>21} | " + "  ".join(f"{v:.6f}" for v in pi))
a, b, c = metodi.values()
print("the three results agree to 1e-10:",
      bool(np.allclose(a, b, atol=1e-10) and np.allclose(b, c, atol=1e-10)))
print("π P = π :", bool(np.allclose(a @ P, a, atol=1e-12)))
print(f"mean return time to / : 1/π₁ = {1 / a[0]:.2f} steps")

Output obtained:

eigenvalues of Pᵀ (complex): +1.0000+0.0000i  +0.7547+0.0000i  -0.6667+0.0000i  -0.2940+0.1549i  -0.2940-0.1549i  -0.5000+0.0000i
chosen index: 0 | eigenvalue: +1.000000+0.000000i
raw eigenvector (real part): +0.5547  +0.5547  +0.3698  +0.3698  +0.1849  +0.2774
        linear system | 0.240000  0.240000  0.160000  0.160000  0.080000  0.120000
     left eigenvector | 0.240000  0.240000  0.160000  0.160000  0.080000  0.120000
         power method | 0.240000  0.240000  0.160000  0.160000  0.080000  0.120000
the three results agree to 1e-10: True
π P = π : True
mean return time to / : 1/π₁ = 4.17 steps

What the code does. The first method builds Pᵀ − I and overwrites its last row with all ones: this is legitimate because the n equations of π = πP are linearly dependent (their sum is 0 = 0), so one of them can be thrown away without losing information. The second calls numpy.linalg.eig on Pᵀ, and here there are three traps in three lines. First: eig on a non-symmetric matrix returns complex numbers, always, even when the eigenvalues are real — in the output four of the six have imaginary part 0.0000, and two are a genuine conjugate pair, −0.2940 ± 0.1549i. Second: eigenvalue 1 is not «the first» in the list, because the ordering is decided by LAPACK and is not guaranteed; you select it with argmin(|λ − 1|), and only then, on that one, do you drop with np.real an imaginary part that is exactly zero.

Third: the eigenvector must be renormalised to sum 1, because NumPy returns it with Euclidean norm 1 — a legitimate vector, but not a probability distribution. The third method starts from the uniform and multiplies a thousand times. The three agree to 1e-10, and the stationary distribution is exactly (6, 6, 4, 4, 2, 3)/25: rational numbers, because the starting probabilities were rational too.

What it tells us. The pages / and /blog share almost half the traffic, and /strumenti gets 8% of it even though the home page links to it exactly as it links to the other two: the difference is that from /strumenti you leave at once and never come back, whereas in the /blog/a1-/blog/a2 pair you bounce. The last line is Kac's theorem: the mean return time to a state is 1/πᵢ, so the crawler comes back to the home page every 4.17 clicks on average. It is also the right way to read the stationary distribution: not «the final probability», but the fraction of time the process spends on each page over a long horizon.

The spectral gap

Knowing that the chain converges does not say when. The answer lies in the second eigenvalue. Order the eigenvalues by decreasing modulus, 1 = |λ₁| > |λ₂| ≥ |λ₃| ≥ …, and assume for simplicity that P is diagonalisable. Then every initial distribution can be written as π₀ = π + Σ_{k≥2} cₖ wₖ, where the wₖ are the left eigenvectors of the eigenvalues λₖ. Applying P once multiplies each piece by its eigenvalue, and applying it n times

πₙ = π + Σ_{k≥2} cₖ λₖⁿ wₖ.

The term π stays put because its eigenvalue is 1; all the others contract, and the last one to die is the largest in modulus. So ‖πₙ − π‖ ≈ C · |λ₂|ⁿ: the error decays geometrically with ratio |λ₂|. The quantity 1 − |λ₂| is called the spectral gap, and the number of steps needed for the error to shrink by a factor of e is about

t_mix ≈ 1 / (1 − |λ₂|).

To drop below a threshold ε you need instead about ln(1/ε) / (1 − |λ₂|) steps: the gap gives the scale, the logarithm of the precision gives the multiplying factor. Let us look at the numbers of our chain.

moduli = sorted(np.abs(np.linalg.eigvals(P)), reverse=True)
print("moduli of the eigenvalues:", "  ".join(f"{m:.4f}" for m in moduli))
lam2 = moduli[1]
print(f"|λ₂| = {lam2:.4f} | gap = 1 − |λ₂| = {1 - lam2:.4f} | t_mix ≈ {1 / (1 - lam2):.2f} steps")

Output obtained:

moduli of the eigenvalues: 1.0000  0.7547  0.6667  0.5000  0.3323  0.3323
|λ₂| = 0.7547 | gap = 1 − |λ₂| = 0.2453 | t_mix ≈ 4.08 steps

What it tells us. Only one eigenvalue has modulus 1 (consistent with Perron-Frobenius: the chain is irreducible and aperiodic), the last two are a complex conjugate pair — perfectly legitimate, because P is not symmetric — and the second in modulus is 0.7547. The gap is wide, the mixing time is four steps: consistent with what we had already seen, namely that after four clicks the distribution had almost forgotten the starting point.

Measured: two sites with the same stationary distribution

The gap cannot be read off the stationary distribution: two chains can have the same limiting distribution and take wildly different times to get there. Let us build two six-page mini-sites, both perfectly symmetric and therefore both with uniform stationary distribution at 1/6. In the first the pages are split into two sections of three, and the visitor stays inside their section 99% of the time. In the second every page links to all the others, like a global menu.

def sito_a_sezioni(interno):
    """Six pages in two sections of three: with probability 'interno' you stay in the section."""
    M = np.zeros((6, 6))
    for i in range(6):
        for j in range(6):
            if i == j:
                continue
            if j // 3 == i // 3:
                M[i, j] = interno / 2.0          # the two sibling pages
            else:
                M[i, j] = (1.0 - interno) / 3.0  # the three pages of the other section
    return M

lento = sito_a_sezioni(0.99)
veloce = np.full((6, 6), 1 / 5) - np.eye(6) / 5  # every page links to all the others

for nome, M in (("closed sections", lento), ("global menu", veloce)):
    moduli = sorted(np.abs(np.linalg.eigvals(M)), reverse=True)
    pi = stazionaria_sistema_lineare(M)
    print(f"{nome:>15} | π₁ = {pi[0]:.4f}, uniform: {bool(np.allclose(pi, 1 / 6))} | "
          f"|λ₂| = {moduli[1]:.4f} | t_mix ≈ {1 / (1 - moduli[1]):.2f}")

partenza = np.zeros(6)
partenza[0] = 1.0
print("   n   closed sections    global menu")
for n in (1, 5, 10, 50, 100):
    da_lento = np.abs(distribuzione_dopo(lento, partenza, n) - 1 / 6).sum()
    da_veloce = np.abs(distribuzione_dopo(veloce, partenza, n) - 1 / 6).sum()
    print(f"{n:>4}   {da_lento:>14.4f}   {da_veloce:>12.4f}")

for nome, M in (("closed sections", lento), ("global menu", veloce)):
    pi_n, passi = partenza.copy(), 0
    while np.abs(pi_n - 1 / 6).sum() >= 0.01:
        pi_n = pi_n @ M
        passi += 1
    print(f"{nome:>15} | steps to fall below 0.01: {passi}")

Output obtained:

closed sections | π₁ = 0.1667, uniform: True | |λ₂| = 0.9800 | t_mix ≈ 50.00
    global menu | π₁ = 0.1667, uniform: True | |λ₂| = 0.2000 | t_mix ≈ 1.25
   n   closed sections    global menu
   1           1.3133         0.3333
   5           0.9039         0.0005
  10           0.8171         0.0000
  50           0.3642         0.0000
 100           0.1326         0.0000
closed sections | steps to fall below 0.01: 228
    global menu | steps to fall below 0.01: 4

What the code does. sito_a_sezioni(0.99) builds the «two almost separate sections» chain: from every page, 99% of the probability is split between the two siblings in the section and 1% among the three pages of the other section. veloce is the complete graph with uniform weights. Both have the same uniform stationary distribution — np.allclose(pi, 1/6) checks it — but their second eigenvalues are 0.98 and 0.20. The table measures the L¹ distance between the distribution at step n and the uniform one, starting from a single page: for the site with closed sections the distance decays exactly like 0.98ⁿ, and after 100 steps it is still at 0.13.

What it tells us. The ratio between the two measured mixing times is 228/4 = 57, whereas the ratio between the 1/(1 − |λ₂|) estimates is 50/1.25 = 40: same order of magnitude, not the same number. The discrepancy is not noise, it is the ln(1/ε) factor the crude estimate ignores: with ε = 0.01 it is 4.6, and 4.6 × 50 = 230, which is practically the 228 we measured. Trade-off. The spectral gap is the cheapest estimate of convergence time — a single eigenvalue — but it is an asymptotic estimate: it describes the behaviour for large n, and it can be far too optimistic in the first steps if the constant C in front is large. When the number of states is large, computing λ₂ costs more than simulating: you use scipy.sparse.linalg.eigs with k=2 and not numpy.linalg.eig, which would require the dense matrix.

Detailed balance and reversibility

There is a condition stronger than stationarity, and far easier to check. A distribution π satisfies detailed balance with respect to P if, for every pair of states,

πᵢ Pᵢⱼ = πⱼ Pⱼᵢ.

In words: the flow of probability going from i to j exactly balances the one coming back, pair by pair. Stationarity asks for much less: that the total inflow and outflow of each state balance, not the individual exchanges. A chain that satisfies detailed balance is called reversible: watched backwards it has the same law as when watched forwards.

Detailed balance implies stationarity. Sum the relation over all possible i: Σᵢ πᵢPᵢⱼ = Σᵢ πⱼPⱼᵢ = πⱼ Σᵢ Pⱼᵢ = πⱼ · 1 = πⱼ. The left-hand side is exactly the j-th component of πP, so πP = π. ∎ The converse is false, and our crawl graph proves it: it has a stationary distribution, but the flow / → /strumenti is 0.24 × 1/3 = 0.08 while the return flow /strumenti → / is 0.08 × 1/2 = 0.04. The books balance only around the full loop, through /contatti.

def scarto_bilancio_dettagliato(P, pi):
    """max |πᵢPᵢⱼ − πⱼPⱼᵢ| over all pairs: zero if and only if the chain is reversible."""
    F = pi[:, None] * P                     # F[i, j] = probability flow from i to j
    return np.abs(F - F.T).max()

pi = stazionaria_sistema_lineare(P)
scarto = scarto_bilancio_dettagliato(P, pi)
print(f"crawl graph: maximum gap = {scarto:.6f} | reversible: {bool(scarto < 1e-12)}")

# same site with reciprocal links: a walk on an UNDIRECTED graph
archi = [(0, 1), (0, 4), (0, 5), (1, 2), (1, 3), (2, 3), (4, 5)]
A = np.zeros((6, 6))
for i, j in archi:
    A[i, j] = A[j, i] = 1.0
gradi = A.sum(axis=1)
R = A / gradi[:, None]
pi_r = gradi / gradi.sum()
print("degrees:", [int(g) for g in gradi], "| π = degree/2m:", "  ".join(f"{v:.6f}" for v in pi_r))
print("π R = π :", bool(np.allclose(pi_r @ R, pi_r)))
print("detailed balance verified to 1e-12:", bool(scarto_bilancio_dettagliato(R, pi_r) < 1e-12))

Output obtained:

crawl graph: maximum gap = 0.040000 | reversible: False
degrees: [3, 3, 2, 2, 2, 2] | π = degree/2m: 0.214286  0.214286  0.142857  0.142857  0.142857  0.142857
π R = π : True
detailed balance verified to 1e-12: True

What it tells us. If the links were all reciprocal — that is, if the graph were undirected — the chain would become reversible and its stationary distribution could be read off without solving any system: πᵢ = dᵢ / Σⱼ dⱼ, the degree of the node divided by the sum of the degrees. The check is immediate: πᵢ Pᵢⱼ = (dᵢ/2m)(1/dᵢ) = 1/2m, which is symmetric in i and j. A closed formula instead of a solve: that is what reversibility buys you.

Why it really matters. Detailed balance is used the other way round. In the lesson on sampling and Monte Carlo methods the problem was drawing samples from a distribution you can only evaluate up to a constant. The solution is to build a Markov chain that has that distribution as its stationary one, run it, and take the visited states as samples. Designing a P that satisfies π = πP directly would be hard; imposing detailed balance, which is a local condition pair by pair, is easy — and that is exactly what the Metropolis-Hastings acceptance ratio does. A chain can be designed to have the stationary distribution you want, and this equation is the tool you design it with.

The experiment: a real surfer, the dangling page and teleportation

Checking by simulation: a million clicks

Everything we have computed so far is linear algebra: π is the solution of a system, and no surfer has ever moved. Yet the stationary distribution has a precise operational meaning — it is the fraction of time that a single surfer, left walking long enough, spends on each page. That is a checkable promise, and we check it: one surfer, a million clicks, fixed seed.

from bisect import bisect_right

def simula(P, partenza, passi, seme=42):
    """Follows a single surfer and counts how many times it visits each page."""
    rng = np.random.default_rng(seme)
    cumulate = [list(riga) for riga in P.cumsum(axis=1)]   # cumulative sums, row by row
    for riga in cumulate:
        riga[-1] = 1.0                                     # 1/3+1/3+1/3 gives 0.9999999999999999
    sorteggi = rng.random(passi)                           # all the random numbers in one go
    conteggi = [0] * len(P)
    stato = partenza
    for u in sorteggi:
        stato = bisect_right(cumulate[stato], u)           # skips zero-probability states by itself
        conteggi[stato] += 1
    return np.array(conteggi) / passi

frequenze = simula(P, pagine.index("/contatti"), 1_000_000)
pi = stazionaria_sistema_lineare(P)
print("       page   frequency    stationary        gap")
for nome, f, p in zip(pagine, frequenze, pi):
    print(f"{nome:>11}   {f:>9.5f}   {p:>11.5f}   {f - p:>+8.5f}")
print(f"maximum gap over 10⁶ steps: {np.abs(frequenze - pi).max():.5f}")

Output obtained (seed 42, NumPy 2.3.5):

       page   frequency    stationary        gap
          /     0.24037       0.24000   +0.00037
      /blog     0.23992       0.24000   -0.00008
   /blog/a1     0.15967       0.16000   -0.00033
   /blog/a2     0.15957       0.16000   -0.00043
 /strumenti     0.08018       0.08000   +0.00018
  /contatti     0.12029       0.12000   +0.00029
maximum gap over 10⁶ steps: 0.00043

What the code does. The trick for not taking half an hour is all in the first few lines: the cumulative sums of each row are computed once, the random numbers are drawn all at once, and choosing the next state becomes a binary search over six numbers. bisect_right also does one thing for free: on zero-probability states the cumulative sum does not grow, so they can never be chosen, and no special case is needed. The seed is fixed at 42, so the simulation is reproducible digit by digit.

What it tells us. The measured visit frequency and the computed stationary distribution agree to within a maximum discrepancy of 0.00043. That is the right order of magnitude: the standard deviation of the empirical frequency over N steps scales like 1/√N, and 1/√(10⁶) = 0.001. Gaining one decimal digit would take a hundred million steps — the same law of diminishing returns that governs every Monte Carlo estimate, and that you saw in the lesson on sampling. Trade-off. On six states simulation is a luxury: the linear system gives the exact answer in microseconds. On millions of states everything flips, because the solve becomes impractical and simulation stays linear in the number of steps: that is why in module 08 returns are estimated by walking, not by solving.

The dangling page problem

So far the mini-site was healthy. Let us break it in the way most common in a real crawl:

link_pozzo = dict(link)
link_pozzo["/contatti"] = ["/contatti"]       # no outgoing link: the crawler stays there
P_pozzo = matrice_di_transizione(pagine, link_pozzo)

descrivi("with the sink", P_pozzo)
pi_pozzo = stazionaria_potenze(P_pozzo, 2000)
print("stationary:", "  ".join(f"{v:.4f}" for v in pi_pozzo))
moduli = sorted(np.abs(np.linalg.eigvals(P_pozzo)), reverse=True)
print(f"|λ₂| = {moduli[1]:.4f} | t_mix ≈ {1 / (1 - moduli[1]):.2f} steps")

Output obtained:

with the sink: 2 class(es), irreducible = False
     transient | ['/', '/blog', '/blog/a1', '/blog/a2', '/strumenti']
     recurrent | ['/contatti']
stationary: 0.0000  0.0000  0.0000  0.0000  0.0000  1.0000
|λ₂| = 0.9278 | t_mix ≈ 13.84 steps

What it tells us. The chain has split into two classes: the five normal pages have become transient (the crawler passes through them finitely many times, then abandons them forever) and /contatti is the only closed class. The stationary distribution still exists and is still unique, but it is degenerate: all the mass on one page, zero on all the others. As a measure of importance it is worth nothing — and note that the spectral gap does not flag the problem, quite the opposite: with |λ₂| = 0.9278 the chain converges slowly, produces numbers that look like they are moving, and still arrives at a useless answer. A small gap is not a bug; a degenerate closed class is.

The classic fix is the teleportation factor of Page, Brin, Motwani and Winograd (1999): with probability d the surfer follows a link, with probability 1 − d they get bored and jump to a random page of the site. The new matrix is

G = d · P + (1 − d) · (1/n) · J,

where J is the all-ones matrix and n the number of pages. Every cell of G is strictly positive, so the chain is irreducible and aperiodic by construction, and Perron-Frobenius guarantees a unique, positive stationary distribution. On top of that there is an exact result on the eigenvalues (Langville and Meyer, 2004): if P has eigenvalues 1, λ₂, λ₃, …, then G has eigenvalues 1, dλ₂, dλ₃, …. Teleportation does not move eigenvalue 1 and multiplies all the others by d: the spectral gap always widens.

def con_teletrasporto(P, d=0.85):
    """G = d·P + (1 − d)/n · J : with probability 1 − d you jump to a random page."""
    n = len(P)
    return d * P + (1 - d) / n * np.ones((n, n))

G = con_teletrasporto(P_pozzo)
descrivi("with teleportation", G)
pi_G = stazionaria_sistema_lineare(G)
for nome, v in zip(pagine, pi_G):
    print(f"{nome:>11} | {v:.4f}")
prima = sorted(np.abs(np.linalg.eigvals(P_pozzo)), reverse=True)[1]
dopo = sorted(np.abs(np.linalg.eigvals(G)), reverse=True)[1]
print(f"|λ₂| before = {prima:.4f} | after = {dopo:.4f} | ratio = {dopo / prima:.4f}")

Output obtained:

with teleportation: 1 class(es), irreducible = True
     recurrent | ['/', '/blog', '/blog/a1', '/blog/a2', '/strumenti', '/contatti']
          / | 0.0888
      /blog | 0.1499
   /blog/a1 | 0.1173
   /blog/a2 | 0.1173
 /strumenti | 0.0502
  /contatti | 0.4765
|λ₂| before = 0.9278 | after = 0.7886 | ratio = 0.8500

What the code does. con_teletrasporto mixes the original matrix with the uniform one, with weight d = 0.85 — the value used in the original PageRank report. The chain becomes irreducible again, the stationary distribution becomes informative again, and the ratio between the second eigenvalues is 0.8500, that is exactly d: the Langville and Meyer result, verified numerically to four digits. The mixing time drops from 13.84 to 4.73 steps, and this explains why PageRank is computed with the power method: with d = 0.85 a few dozen iterations are enough, on any graph.

The limit of the fix, declared. Look at the number for /contatti: 0.4765, almost half of all the mass. Teleportation has repaired the mathematics — existence, uniqueness, convergence — but it has not taken away the dangling page's structural advantage: you enter it from two sides and never leave except out of boredom. The complete remedy, in the 1999 report, is to handle pages with no outgoing links by redistributing their mass over all pages, instead of letting them cycle on themselves. The rule to take away. Teleportation is a change to the model, not a numerical patch: it changes the chain whose stationary distribution you are computing, and whoever reads the result has the right to know that d = 0.85 is a choice of yours.

From process to decision

So far the process was something you underwent. The crawler does not choose: the matrix P is given, and it follows it. Write the transition as P(s′ | s) — the probability of ending in state s′ starting from s — and you have described everything that can happen.

Now add one single thing: in every state you can take an action, chosen from a finite set, and the action changes the transition matrix. The transition becomes

P(s′ | s, a),

that is, no longer one matrix but one matrix per action. This object is called a Markov decision process, and a rule assigning an action to every state is called a policy: once a policy is fixed, the matrices collapse back into one, and you are in a Markov chain again, with its stationary distribution and its spectral gap. That is why the chain comes first: a decision process is a family of chains, indexed by policies, plus a criterion for choosing one.

From here on module 08 adds two ingredients that are not needed here: a reward attached to the transitions and a criterion for comparing policies over a long horizon. The Bellman equation is born from there. But the object all of this rests on — the state space, the probabilistic transition, the property that the future depends only on the present — is the one you built in this lesson, and the vocabulary is the same: state, transition, distribution, convergence.

Typical mistakes (symptom, cause, remedy)

1. Row/column convention swapped. Symptom: π does not add up, or it «converges» to a distribution that has nothing to do with the graph. Cause: you wrote P @ pi with P row-stochastic, that is, you simulated the transposed chain. Remedy: declare the convention in a comment, check P.sum(axis=1) and not P.sum(axis=0), and verify the first step by hand on a state with a single outgoing link.
2. Looking for a convergence that does not exist. Symptom: the power method oscillates between two or more values and never settles. Cause: the chain is periodic, and Perron-Frobenius guarantees convergence only with aperiodicity. Remedy: compute the period, or look at how many eigenvalues have modulus 1; if there is more than one, the stationary distribution exists but is not a limit, and must be computed by solving the linear system.
3. Eigenvector not normalised to sum 1. Symptom: the components of π add up to 1.7 or to −1. Cause: numpy.linalg.eig returns eigenvectors of Euclidean norm 1, with the sign decided by the algorithm. Remedy: divide by the sum, not by the norm, and if the sum is negative flip the sign of the whole vector before dividing.
4. Right eigenvector instead of the left one. Symptom: you get a vector with all equal components, which looks like a plausible uniform distribution. Cause: you called eig(P) instead of eig(P.T), and the right eigenvector of eigenvalue 1 is always the all-ones vector. Remedy: if the result is uniform on a chain that is not symmetric, you have picked the wrong side.
5. Treating a non-Markovian state as Markovian. Symptom: the model predicts badly, and worst where history matters (queues, delays, penalties). Cause: the state is too poor. Remedy: list the pieces of information you would use to make the prediction by hand; if you use one that is not in the state, the state has to be widened — and declare how much the space grows.
6. Power method stopped too early. Symptom: two runs with a different number of iterations give different answers at the third digit. Cause: the spectral gap is small, and after n steps the error is still |λ₂|ⁿ. Remedy: do not fix the iterations, fix a threshold on the change between two consecutive steps, and as a check compute |λ₂| once to know what order of magnitude to expect.
7. Reading the stationary distribution as «the final answer». Symptom: one concludes the system «will end up» in that state. Cause: confusion between limiting distribution and destination. Remedy: the stationary distribution is the fraction of time spent in each state over the long run — the process keeps moving, always.

The seven mistakes are recognised by their symptom: four give wrong but plausible numbers (convention, right eigenvector, normalisation, powers stopped early) and three are model errors (periodicity ignored, state too poor, stationary distribution read as a destination). Against the first four it is enough to compute π in two ways and compare; against the three model errors you need to look at the graph before the formulas.

Best practice

• Declare the convention — row or column — in the first comment line, and check it with an assert on the sums: it costs one line and prevents the most expensive mistake.
• Check irreducibility and period before looking for a stationary distribution, not after: they are two graph traversals and a gcd, and they tell you whether the question you are asking has an answer.
• Compute π in two independent ways and compare: the linear system and the power method use different algorithms, and a disagreement is almost always an error in building the matrix.
• Look at the second eigenvalue before choosing the number of iterations, and use a threshold on the change instead of a fixed number of steps.
• Declare the model's parameters — the teleportation factor, the uniformity assumption on links, time homogeneity — next to the result: they are choices of yours, not properties of the data.
• For large sparse matrices use scipy.sparse.linalg.eigs with a small k, and never build Pⁿ: multiply one vector at a time.

The six rules are checks of very different cost: the convention and the two ways of computing cost one line each, the check of irreducibility and period two visits of the graph, the declaration of parameters is one sentence. If you must pick only one, pick the comparison between two independent methods: it exposes almost all errors in building the matrix.

Minimal glossary

• Markov property: the future depends on the past only through the present state.
• Transition matrix: Pᵢⱼ = probability of going from i to j; non-negative rows that sum to 1.
• Irreducible: from every state you can reach every other state; a single communication class.
• Aperiodic: the greatest common divisor of the return times is 1; without this property the distribution oscillates.
• Recurrent / transient: you come back with probability 1 / you come back only finitely many times.
• Stationary distribution: π = πP; left eigenvector of eigenvalue 1, normalised to sum 1.
• Spectral gap: 1 − |λ₂|; it governs the speed at which the error |λ₂|ⁿ dies out.
• Mixing time: number of steps needed to come within a threshold of the stationary distribution; about ln(1/ε)/(1 − |λ₂|).
• Detailed balance: πᵢPᵢⱼ = πⱼPⱼᵢ; it implies stationarity, not the other way round.
• Kac's theorem: the mean return time to a state is 1/πᵢ.

The ten entries follow the order of the sections and fall into three groups: the chain (property, matrix), the structure of the graph (irreducible, aperiodic, recurrent and transient) and the dynamics (stationary, spectral gap, mixing, detailed balance, Kac). After the definition of stationary distribution the other entries are its consequences or its measuring tools.

Mini-project

Build a module that, given a dictionary of links, answers every question of this lesson.

1. Write analizza(link) returning: the transition matrix, the communication classes with their recurrent/transient label, the period, the stationary distribution computed in two independent ways and |λ₂|;
2. apply it to three graphs: the one in the lesson, the same with /contatti having no outgoing links, and the same with teleportation at d = 0.85; print the five answers for each in a table;
3. handle and document the two edge cases: a page that links only to itself (a closed class of a single state) and a graph with two separate components (stationary distribution not unique: your code has to notice and say so, instead of returning a random number).

The three requests build a single module that must always answer the same five questions on any graph: the first defines the interface, the second tests it on three cases, the third pushes it to the edge cases where the Perron-Frobenius hypotheses fail. The edge cases are what distinguishes a script from a module.

Deliverable: a script that prints the table for the three chains and, at the bottom, a line declaring which assumption of the model you consider the least defensible, and why.

Challenge

• Easy. Pen and paper: given P = [[0.9, 0.1], [0.5, 0.5]], compute π by solving π = πP with π₁ + π₂ = 1. Solution: the first equation is π₁ = 0.9π₁ + 0.5π₂, that is 0.1π₁ = 0.5π₂, that is π₁ = 5π₂. With the constraint 5π₂ + π₂ = 1 you get π₂ = 1/6 and π₁ = 5/6. Check with the second equation: 0.1·(5/6) + 0.5·(1/6) = 1/12 + 1/12 = 1/6 ✓.
• Medium. Pen and paper: prove that detailed balance implies stationarity, then build a counterexample to the converse with three states (hint: a directed cycle 1 → 2 → 3 → 1 with all probabilities equal to 1 has a uniform stationary distribution, but the flows in one direction do not balance those in the other).
• Hard. Code: build two four-state chains with the same uniform stationary distribution, one periodic of period 2 and one aperiodic, and show the difference numerically: for the first the power method oscillates, for the second it converges; count the eigenvalues of modulus 1 in both cases.
• Expert. Code: estimate the mixing time by simulation — start from every possible state, measure after how many steps the L¹ distance from the stationary distribution drops below 0.01, take the maximum over starting states — and compare it with 1/(1 − |λ₂|) and with ln(100)/(1 − |λ₂|) on at least five different chains. Which of the two estimates is systematically better?
• Enterprise. Take a state of a system from your own work that evolves in discrete steps — the states of a ticket, of an order, of a lead in the funnel — and estimate P from historical data by counting transitions (P̂ᵢⱼ = nᵢⱼ / Σₖ nᵢₖ). Report the stationary distribution, the mean time spent in each state and |λ₂|. Then the hard part: check whether the Markov property really holds, by comparing the distribution of the transitions leaving a state as a function of the previous state. If the two distributions differ, the state has to be widened, and your one-step model is lying.

The five levels climb from the computation by hand to an analysis on your own data: the first two are solved with pen and paper (the stationary distribution of a 2×2 matrix and the link between detailed balance and stationarity), the third and the fourth ask for code to compare periodic and aperiodic chains and mixing-time estimates, the last estimates P from your data and checks that the Markov property really holds.

The code exercises of this lesson are in a single file, esercizi.py (about 11 KB). Four code exercises on the themes of the lesson: building the transition matrix from a dictionary of links, computing the stationary distribution by solving the linear system, verifying it by simulating a surfer for a million steps, and estimating the mixing time by comparing it with the spectral gap. A skeleton with the parts to complete and, at the end, the commented solutions. It requires Python 3.12+ and NumPy 2.3. Run it with python3 esercizi.py: it works even with the parts to complete left empty and tells you which ones are missing.

Recap

Mental map. Markov property: the future depends on the past only through the state — and you choose the state. P row-stochastic, distribution as a row vector, πₙ = π₀Pⁿ. Irreducible ⇒ the stationary distribution exists, is unique, is positive. Aperiodic ⇒ you get there. π = πP ⟺ Pᵀπᵀ = πᵀ: left eigenvector of eigenvalue 1. The error decays like |λ₂|ⁿ, the mixing time is about ln(1/ε)/(1 − |λ₂|). Detailed balance πᵢPᵢⱼ = πⱼPⱼᵢ ⇒ stationarity, not the other way round. A chain is something you undergo; with P(s′ | s, a) you decide.

Cheat sheet.

3 concepts: the state summarises the past, it does not erase it; the stationary distribution is the left eigenvector of eigenvalue 1; the second eigenvalue governs the speed. 3 mistakes: row/column convention swapped; looking for a limit on a periodic chain; using the right eigenvector instead of the left one. 3 best practices: check irreducibility and period before computing π; compute it in two independent ways; stop the power method on a threshold, not on a fixed number of steps. 1 exercise: the mini-project of this lesson.

Next lesson: the mathematics ends here. From the next lesson on, the models: the process of a machine learning project, and the three metrics that lie when the classes are imbalanced.

Questions to take further

• How do you estimate a transition matrix from data when some transitions have never been observed, and how do you avoid declaring them impossible? (Laplace correction and Dirichlet priors: lesson 13 of this module.)
• What happens when the state space is continuous or infinite, and the matrix becomes a transition kernel? (That is the case of Markov chain Monte Carlo methods: lesson 14 of this module.)
• How do you choose the best action when each state offers many and the result only shows up after many steps? (Markov decision process and the Bellman equation, module 08.)

The three questions are hooks to other lessons and each already has an answer in brackets: the first concerns transitions never observed and points to the Laplace correction, the second continuous state spaces and Markov chain Monte Carlo methods, the third the choice of actions in the Markov decision process. They help you decide what to read next.