30 - Catene di Markov: la proprietà che dimentica il passato, e la distribuzione che non cambia più

Corso Machine Learning Engineer · Modulo 02 · Matematica e statistica · Lezione 15

Obiettivo. Riconoscere quando un processo è una catena di Markov, scriverne la matrice di transizione e calcolarne la distribuzione stazionaria in tre modi indipendenti. Alla fine sai dire se una catena ha una distribuzione limite, sai calcolarla come autovettore sinistro di autovalore 1, sai stimare quanti passi servono perché si stabilizzi guardando il secondo autovalore, e sai riconoscere la differenza fra un processo che subisci e un processo in cui decidi.

Falso mito da sfatare. «Markov vuol dire che il passato non conta.» No: la proprietà di Markov dice che il passato non conta in più rispetto al presente, perché lo stato — se è scelto bene — lo riassume già tutto. Non è una proprietà del fenomeno: è una proprietà della coppia fenomeno-stato. Scegliere male lo stato è il modo più comune di costruire un modello markoviano sbagliato.

Introduzione

Questa è l'ultima lezione del modulo di matematica, e chiude due cerchi insieme. Il primo: nella lezione sugli autovalori, la SVD e la PCA hai imparato che certe matrici hanno direzioni speciali, e lì gli autovettori servivano a comprimere i dati. Qui servono a tutt'altro: l'autovettore di autovalore 1 di una matrice di transizione è la distribuzione a cui un processo tende, e il secondo autovalore dice quanto ci mette ad arrivarci. È l'unico uso serio degli autovalori, in questo corso, fuori dalla PCA.

Il secondo cerchio guarda avanti. Nel modulo 08 il processo decisionale di Markov è il linguaggio in cui si scrive ogni problema di apprendimento per rinforzo, e l'equazione di Bellman ne è la conseguenza. Nessuna delle due lezioni, però, definisce la catena che sta sotto: la matrice di transizione, la stazionaria, la condizione perché esista. La costruiamo qui, così che là resti solo da aggiungere la decisione.

Risultato finale di questa lezione. Dato un grafo di sei pagine con i link dichiarati, sai costruire la matrice di transizione, calcolare la distribuzione stazionaria in tre modi che coincidono a 1e-10, stimare il tempo di stabilizzazione dal gap spettrale, diagnosticare una pagina senza link uscenti e correggerla con il fattore di teletrasporto.

Roadmap della lezione.

1. La proprietà di Markov, e i processi che non la rispettano finché non allarghi lo stato.
2. La matrice di transizione, la convenzione riga-colonna e il grafo di crawl a sei pagine.
3. La distribuzione dopo n passi, Chapman-Kolmogorov e il conto di Pⁿ per quadrati ripetuti.
4. Classi di comunicazione, stati ricorrenti e transitori.
5. Irriducibilità, periodo e aperiodicità, con la catena che oscilla per sempre.
6. La distribuzione stazionaria come autovettore sinistro, Perron-Frobenius, tre modi di calcolarla.
7. Il gap spettrale, la velocità di convergenza e il tempo di mescolamento, misurati.
8. Il bilancio dettagliato, e perché su di esso si progettano le catene.
9. Il problema del pozzo e il fattore di teletrasporto.
10. Dalla catena al processo decisionale: la riga che porta al modulo 08.

La scaletta ha due metà: le prime cinque voci costruiscono il vocabolario (proprietà, matrice, potenze, classi, periodo), le ultime cinque lo usano per rispondere alle domande pratiche (qual è la stazionaria, quanto ci mette, come si progetta una catena, che cosa succede se si rompe). Se ti interessa solo l'apprendimento per rinforzo, la lezione si può leggere fino alla stazionaria e poi saltare all'ultima voce.

La proprietà di Markov

L'analogia: il gioco dell'oca

Sei alla casella 34 del gioco dell'oca. Per sapere dove finirai al prossimo lancio ti basta una cosa sola: che sei alla 34. Non serve sapere se ci sei arrivato con un 6 e un 2, con tre lanci sfortunati o con una scorciatoia: il tabellone non se lo ricorda, e il dado nemmeno. Il futuro dipende dal passato solo attraverso la posizione attuale.

Un processo fatto così si chiama catena di Markov, dal matematico russo Andrej Andreevič Markov, che nel 1906 li introdusse per estendere la legge dei grandi numeri a variabili dipendenti fra loro. È la struttura più semplice che non sia l'indipendenza: non «ogni passo è scollegato dagli altri», ma «ogni passo è collegato solo al precedente».

La definizione formale

Definizione. Una successione di variabili aleatorie X₀, X₁, X₂, … a valori in un insieme finito di stati è una catena di Markov se, per ogni n e per ogni storia,

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

Leggila come una frase: la probabilità di finire in j al passo successivo, sapendo tutta la storia, è la stessa che si ha sapendo soltanto dove si è adesso. Se questa probabilità non dipende neanche da n — cioè le regole del gioco non cambiano nel tempo — la catena si dice omogenea, e scriviamo Pᵢⱼ = P(Xₙ₊₁ = j | Xₙ = i). Tutte le catene di questa lezione sono omogenee, ed è bene dichiararlo: è un'ipotesi, non un dato di natura.

L'avvertenza che evita l'errore classico. La proprietà non dice che il passato è irrilevante. Dice che il passato è già riassunto nello stato. Sono due affermazioni molto diverse: la prima è quasi sempre falsa nei sistemi reali, la seconda è una richiesta di progettazione che tocca a te soddisfare, scegliendo cosa mettere dentro lo stato.

Due processi che non sono markoviani (finché non allarghi lo stato)

Primo esempio: la navigazione con il pulsante «indietro». Un visitatore si muove su un sito e ogni tanto preme «indietro», che lo riporta alla pagina da cui era venuto. Se prendi come stato «la pagina su cui mi trovo», il processo non è markoviano: la probabilità di finire su /blog dipende da dove stavi prima, e quell'informazione nello stato non c'è. La correzione non è abbandonare Markov, è allargare lo stato: prendi come stato la coppia (pagina corrente, pagina precedente). Sulle coppie la proprietà vale di nuovo, al prezzo di passare da 6 stati a (al più) 36.

Secondo esempio: la somma degli ultimi due lanci di un dado. Sia Sₙ = Dₙ + Dₙ₋₁. Se Sₙ = 7, la distribuzione di Sₙ₊₁ = Dₙ₊₁ + Dₙ dipende da quanto valeva Dₙ: con Dₙ = 6 la somma successiva sta fra 7 e 12, con Dₙ = 1 fra 2 e 7. Lo stato «somma» butta via proprio il pezzo di passato che serve. Correzione: usa come stato Dₙ (o la coppia (Dₙ₋₁, Dₙ) se vuoi continuare a leggere la somma dallo stato).

La regola da portare via. «Non è markoviano» quasi sempre significa «lo stato è troppo povero». Allargare lo stato costa: il numero di stati esplode, e la matrice di transizione cresce con il quadrato. La scelta dello stato è il vero progetto di un modello markoviano, e si paga in memoria.

La matrice di transizione

Con un numero finito di stati, tutte le probabilità di un passo stanno in una matrice quadrata P, dove Pᵢⱼ è la probabilità di passare da i a j. Due proprietà la definiscono: Pᵢⱼ ≥ 0 e, per ogni riga,

Σⱼ Pᵢⱼ = 1.

In parole: da qualunque stato, da qualche parte si finisce. Una matrice con queste due proprietà si dice stocastica per righe. Non è una convenzione universale: metà della letteratura di fisica e di algebra lineare usa matrici stocastiche per colonne, con Pᵢⱼ uguale alla probabilità di andare da j a i.

La convenzione, dichiarata una volta per tutte. In questa lezione P è stocastica per righe, e la distribuzione di probabilità è un vettore riga: se π₀ è la distribuzione iniziale, dopo un passo la distribuzione è π₁ = π₀ P, e dopo n passi

πₙ = π₀ Pⁿ.

Scrivere P @ pi invece di pi @ P è l'errore numero uno di chi implementa una catena per la prima volta: non solleva nessuna eccezione, perché la matrice è quadrata, e restituisce semplicemente i numeri di un'altra catena (quella trasposta). Dichiara la convenzione in un commento sopra la matrice e verificala subito, controllando che le righe sommino a 1 e non le colonne.

Il grafo di crawl a sei pagine

L'esempio portante della lezione è un mini-sito di sei pagine, con i link dichiarati a mano: nessun dataset esterno, nessuna rete scaricata. Un crawler parte da una pagina e a ogni passo segue un link uscente scelto a caso con probabilità uniforme. Le sei pagine sono /, /blog, /blog/a1, /blog/a2, /strumenti, /contatti, e la struttura dei link è quella della figura.

La traduzione in codice è meccanica: un dizionario che per ogni pagina elenca i link uscenti, e una funzione che lo trasforma in matrice distribuendo 1/k sui k link.

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] = probabilità di passare dalla pagina i alla pagina j in un clic."""
    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("ogni riga somma a 1:", bool(np.allclose(P.sum(axis=1), 1.0)))

Output ottenuto (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
ogni riga somma a 1: True

Cosa fa il codice. Il dizionario link è la sola cosa scritta a mano: tutto il resto è derivato. matrice_di_transizione mette 1/len(uscite) nella casella (sorgente, destinazione), che è l'ipotesi «il crawler sceglie un link a caso, tutti con la stessa probabilità» — un'ipotesi, non un dato: un crawler vero segue i link con probabilità diverse, e un visitatore umano ancora di più. La riga di /contatti ha un solo 1.000 perché quella pagina ha un unico link uscente. Il controllo finale sulle righe è l'unica difesa contro un errore di battitura nel dizionario.

La distribuzione dopo n passi

Se la distribuzione di oggi è il vettore riga π₀, quella di domani è π₀P. La riga i di P dice dove va la massa che si trova in i, e il prodotto riga-per-matrice somma i contributi di tutte le partenze. Ripetendo, dopo n passi si ha πₙ = π₀Pⁿ: la potenza n-esima della matrice contiene tutte le probabilità di transizione in n passi, perché (Pⁿ)ᵢⱼ = P(Xₙ = j | X₀ = i).

def distribuzione_dopo(P, pi0, n):
    """Vettore riga per matrice, n volte: π₀ Pⁿ."""
    pi = pi0.astype(float).copy()
    for _ in range(n):
        pi = pi @ P                      # riga per matrice: NON P @ pi
    return pi

pi0 = np.zeros(6)
pi0[pagine.index("/contatti")] = 1.0    # il crawler parte da /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 ottenuto:

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

Cosa fa il codice. Parte dalla certezza di essere su /contatti e applica la matrice una volta per passo. I primi passi sono deterministici o quasi: da /contatti si va per forza su /, da lì si va con un terzo ciascuno su /blog, /strumenti, /contatti. Dopo quattro passi la memoria del punto di partenza è già quasi sparita, e dopo cento la distribuzione non si muove più di 1e-13: è la stazionaria. Attenzione: «non si muove più» non è un fatto automatico, è una conseguenza di due proprietà della catena (irriducibilità e aperiodicità) che tra poco definiremo e che, se mancano, fanno saltare tutto.

Chapman-Kolmogorov, e il conto di Pⁿ

Andare da i a j in m + n passi significa passare per un qualche stato intermedio k dopo m passi. Sommando su tutti i k possibili si ottiene l'equazione di Chapman-Kolmogorov:

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

che è, parola per parola, la definizione di prodotto fra matrici: P^(m+n) = Pᵐ Pⁿ. Non è un teorema profondo, è la proprietà di Markov scritta in notazione matriciale. Ha però una conseguenza pratica: per calcolare P¹⁰⁰ non servono 99 moltiplicazioni, ne bastano una decina se si usa l'esponenziazione per quadrati ripetuti, che spezza l'esponente in potenze di due.

def potenza_veloce(M, n):
    """Mⁿ per quadrati ripetuti; restituisce anche quante moltiplicazioni sono servite."""
    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³ (stesso risultato)  :", bool(np.allclose(P7, P4 @ P3)))
P100, conta = potenza_veloce(P, 100)
print(f"P¹⁰⁰: {conta} moltiplicazioni invece di 99")
print("tutte le righe di P¹⁰⁰ sono uguali:", bool(np.allclose(P100, np.ones((6, 1)) @ P100[:1])))

Output ottenuto:

Chapman-Kolmogorov, P⁷ = P³ P⁴ : True
P⁷ = P⁴ P³ (stesso risultato)  : True
P¹⁰⁰: 9 moltiplicazioni invece di 99
tutte le righe di P¹⁰⁰ sono uguali: True

Che cosa ci dice. L'ultima riga è il fatto più importante delle quattro: in P¹⁰⁰ tutte le righe coincidono, cioè dopo cento passi non conta più da dove sei partito. Una matrice di questo tipo ha rango 1 nel limite, e la riga ripetuta è esattamente la distribuzione stazionaria. Trade-off. I quadrati ripetuti convengono per un n grande su una matrice piccola e densa; se la matrice è enorme e sparsa (un grafo web vero) le potenze la riempiono di zeri che diventano non-zeri, e conviene invece moltiplicare ripetutamente un singolo vettore, che costa un prodotto sparso per passo e non costruisce mai Pⁿ.

Classificazione degli stati

Prima di chiedersi «dove va a finire» una catena bisogna sapere com'è fatto il suo grafo. Quattro definizioni bastano. Lo stato j è accessibile da i (si scrive i → j) se esiste n ≥ 0 con (Pⁿ)ᵢⱼ > 0, cioè se dal grafo si può andare da i a j seguendo le frecce. Due stati comunicano se i → j e j → i. La comunicazione è una relazione di equivalenza, quindi spezza gli stati in classi di comunicazione. Una classe è chiusa se nessuna freccia ne esce.

Da qui seguono le due etichette che contano. Uno stato è ricorrente se, partendo da lì, la probabilità di tornarci prima o poi è 1; è transitorio se quella probabilità è minore di 1, e in quel caso la catena ci torna solo un numero finito di volte, poi lo abbandona per sempre. Su uno spazio di stati finito vale una scorciatoia comoda: gli stati di una classe chiusa sono ricorrenti, quelli di una classe non chiusa sono transitori. La catena è irriducibile se c'è una sola classe, cioè se da ogni pagina si può raggiungere ogni altra pagina.

def raggiungibili(P, i):
    """Stati raggiungibili da i in un numero qualunque di passi (visita del grafo)."""
    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):
    """Una classe è chiusa se nessuna transizione ne esce: su stati finiti, è ricorrente."""
    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)} classe/i, irriducibile = {len(classi) == 1}")
    for classe in classi:
        tipo = "ricorrente" if chiusa(M, classe) else "transitoria"
        print(f"   {tipo:>11} | {[pagine[i] for i in classe]}")

descrivi("grafo di crawl", P)

Output ottenuto:

grafo di crawl: 1 classe/i, irriducibile = True
    ricorrente | ['/', '/blog', '/blog/a1', '/blog/a2', '/strumenti', '/contatti']

Cosa fa il codice. raggiungibili è una visita del grafo che ignora i valori delle probabilità e guarda solo quali sono diverse da zero: per l'accessibilità conta la struttura, non i numeri. classi_di_comunicazione mette insieme gli stati mutuamente raggiungibili; chiusa controlla che non ci siano frecce in uscita. Il nostro grafo ha una sola classe e quindi è irriducibile: ogni pagina è ricorrente, e il crawler passerà infinite volte da ciascuna. Questa stessa funzione descrivi la riuseremo più avanti sulla versione rotta del sito, dove la risposta cambia.

Irriducibilità, periodo e aperiodicità

L'irriducibilità da sola non basta a garantire che la distribuzione si stabilizzi. Serve una seconda proprietà. Il periodo di uno stato i è

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

cioè il massimo comun divisore delle lunghezze dei cammini che tornano in i. Se d(i) = 1 lo stato si dice aperiodico. In una catena irriducibile tutti gli stati hanno lo stesso periodo, quindi si può parlare del periodo della catena. Un auto-anello qualsiasi (Pᵢᵢ > 0) rende la catena aperiodica all'istante, perché mette il numero 1 nell'insieme dei tempi di ritorno; ma non è necessario: bastano due cicli di lunghezze coprime, come il nostro / → /strumenti → / di lunghezza 2 e / → /strumenti → /contatti → / di lunghezza 3.

from math import gcd

def periodo(P, i, massimo=24):
    """gcd delle lunghezze n per cui è possibile tornare in i in n passi."""
    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("periodo di ogni pagina del grafo di crawl:", [periodo(P, i) for i in range(6)])

# due pagine che si linkano solo a vicenda: catena periodica di periodo 2
Q = np.array([[0.0, 1.0],
              [1.0, 0.0]])
print("periodo delle due pagine:", [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) soddisfa π = πQ:", bool(np.allclose(pi_stazionaria @ Q, pi_stazionaria)))
print("moduli degli autovalori:", "  ".join(f"{m:.4f}" for m in sorted(abs(np.linalg.eigvals(Q)), reverse=True)))

Output ottenuto:

periodo di ogni pagina del grafo di crawl: [1, 1, 1, 1, 1, 1]
periodo delle due pagine: [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) soddisfa π = πQ: True
moduli degli autovalori: 1.0000  1.0000

Che cosa ci dice. La catena a due pagine è irriducibile — si va dall'una all'altra e viceversa — eppure la distribuzione non converge mai: oscilla fra (1, 0) e (0, 1) per sempre. Il vettore (0.5, 0.5) è stazionario, nel senso che soddisfa π = πQ, ma nessuna distribuzione iniziale diversa da lui ci arriva. La spia numerica è l'ultima riga: ci sono due autovalori di modulo 1 (sono 1 e −1), e quell'autovalore −1 è esattamente l'oscillazione. La regola da portare via. Irriducibile garantisce che la stazionaria esista e sia unica; aperiodica garantisce che ci si arrivi. Sono due richieste separate, e la seconda si dimentica più spesso.

La distribuzione stazionaria

Definizione. Una distribuzione di probabilità π (vettore riga, componenti non negative che sommano a 1) è stazionaria se

π = π P.

In parole: se oggi la probabilità di trovarsi su ciascuna pagina è π, domani sarà ancora π. Non vuol dire che il crawler sta fermo — continua a saltare da una pagina all'altra — vuol dire che il flusso in entrata di ogni pagina pareggia esattamente quello in uscita. È un equilibrio dinamico, come il livello di una vasca con rubinetto e scarico aperti.

È un autovettore sinistro, e questo chiude il cerchio con gli autovalori

Trasponendo entrambi i membri di π = πP si ottiene

Pᵀ πᵀ = πᵀ,

cioè: πᵀ è un autovettore di Pᵀ relativo all'autovalore 1. Equivalentemente, π è un autovettore sinistro di P con autovalore 1. L'autovalore 1 c'è sempre, per un motivo di una riga: il vettore colonna di soli uni 𝟙 soddisfa P𝟙 = 𝟙 (è la condizione di righe che sommano a 1), quindi 1 è autovalore di P, e una matrice e la sua trasposta hanno gli stessi autovalori. Il vettore 𝟙 è l'autovettore destro; quello sinistro, che non gli assomiglia per niente, è la stazionaria.

Qui si vede perché nella lezione sugli autovalori abbiamo insistito sulla distinzione fra autovettori destri e sinistri. Nella PCA la matrice era simmetrica e la distinzione spariva: per una matrice simmetrica destro e sinistro coincidono. Una matrice di transizione, invece, non è quasi mai simmetrica — /contatti punta a / con probabilità 1, / punta a /contatti con probabilità un terzo — e scambiare i due autovettori è un errore che produce numeri plausibili e sbagliati.

Perron-Frobenius, enunciato con precisione

Teorema (Perron 1907, Frobenius 1912), caso delle matrici stocastiche. Sia P la matrice di transizione di una catena di Markov finita e irriducibile. Allora: (a) esiste una e una sola distribuzione di probabilità π con π = πP, e tutte le sue componenti sono strettamente positive; (b) tutti gli autovalori hanno modulo ≤ 1, e l'autovalore 1 è semplice, cioè con molteplicità 1; (c) se inoltre la catena è aperiodica, allora 1 è l'unico autovalore di modulo 1, e πₙ → π per qualunque distribuzione iniziale π₀.

Lo enunciamo senza dimostrarlo: la dimostrazione standard sta in Norris, Markov Chains, capitolo 1, e la versione generale per matrici non negative in Horn e Johnson, Matrix Analysis, capitolo 8. Lo verificheremo numericamente, che è un'altra cosa e vale meno di una dimostrazione, ma serve a controllare di avere capito l'enunciato. Vale la pena notare dove ciascuna ipotesi entra in gioco: senza irriducibilità la stazionaria può non essere unica (due classi chiuse separate ne danno infinite, tutte le loro combinazioni); senza aperiodicità esiste ed è unica, ma non la si raggiunge, come nell'esempio delle due pagine.

Tre modi di calcolarla

Il primo è il più diretto: π = πP è un sistema lineare omogeneo π(P − I) = 0, che da solo ha infinite soluzioni (tutti i multipli di π). Si rende unico sostituendo una qualunque delle equazioni — che è ridondante — con il vincolo Σπᵢ = 1. Il secondo chiede l'autovettore a una libreria. Il terzo moltiplica un vettore per la matrice finché non si muove più.

def stazionaria_sistema_lineare(P):
    """Risolve (Pᵀ − I)πᵀ = 0 sostituendo l'ultima equazione con Σπᵢ = 1."""
    n = len(P)
    A = P.T - np.eye(n)
    A[-1, :] = 1.0                      # l'equazione ridondante diventa la normalizzazione
    b = np.zeros(n)
    b[-1] = 1.0
    return np.linalg.solve(A, b)

def stazionaria_autovettore(P, mostra=False):
    """Autovettore sinistro di P di autovalore 1 = autovettore destro di Pᵀ."""
    valori, vettori = np.linalg.eig(P.T)                 # attenzione: numeri COMPLESSI
    if mostra:
        print("autovalori di Pᵀ (complessi):",
              "  ".join(f"{v.real:+.4f}{v.imag:+.4f}i" for v in valori))
    k = int(np.argmin(np.abs(valori - 1.0)))             # il più vicino a 1, non "il primo"
    if mostra:
        print(f"indice scelto: {k} | autovalore: {valori[k].real:+.6f}{valori[k].imag:+.6f}i")
        print("autovettore grezzo (parte reale):", "  ".join(f"{v:+.4f}" for v in np.real(vettori[:, k])))
    v = np.real(vettori[:, k])                           # scartare la parte immaginaria ~0
    return v / v.sum()                  # senza questa riga la norma è 1, non la somma

def stazionaria_potenze(P, passi=1000):
    """Metodo delle potenze sul vettore riga, partendo dalla distribuzione uniforme."""
    pi = np.full(len(P), 1.0 / len(P))
    for _ in range(passi):
        pi = pi @ P
    return pi

metodi = {
    "sistema lineare": stazionaria_sistema_lineare(P),
    "autovettore sinistro": stazionaria_autovettore(P, mostra=True),
    "metodo delle potenze": 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("i tre risultati coincidono a 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"tempo medio di ritorno su / : 1/π₁ = {1 / a[0]:.2f} passi")

Output ottenuto:

autovalori di Pᵀ (complessi): +1.0000+0.0000i  +0.7547+0.0000i  -0.6667+0.0000i  -0.2940+0.1549i  -0.2940-0.1549i  -0.5000+0.0000i
indice scelto: 0 | autovalore: +1.000000+0.000000i
autovettore grezzo (parte reale): +0.5547  +0.5547  +0.3698  +0.3698  +0.1849  +0.2774
      sistema lineare | 0.240000  0.240000  0.160000  0.160000  0.080000  0.120000
 autovettore sinistro | 0.240000  0.240000  0.160000  0.160000  0.080000  0.120000
 metodo delle potenze | 0.240000  0.240000  0.160000  0.160000  0.080000  0.120000
i tre risultati coincidono a 1e-10: True
π P = π : True
tempo medio di ritorno su / : 1/π₁ = 4.17 passi

Cosa fa il codice. Il primo metodo costruisce Pᵀ − I e ne sovrascrive l'ultima riga con tutti uno: è lecito perché le n equazioni di π = πP sono linearmente dipendenti (la loro somma è 0 = 0), quindi se ne può buttare una senza perdere informazione. Il secondo chiama numpy.linalg.eig su Pᵀ, e qui ci sono tre trappole in tre righe. Primo: eig su una matrice non simmetrica restituisce numeri complessi, sempre, anche quando gli autovalori sono reali — nell'output quattro dei sei hanno parte immaginaria 0.0000, e due sono una coppia coniugata vera, −0.2940 ± 0.1549i. Secondo: l'autovalore 1 non è «il primo» della lista, perché l'ordine lo decide LAPACK e non è garantito; si seleziona con argmin(|λ − 1|), e solo dopo, su quello, si scarta con np.real una parte immaginaria che vale esattamente zero. Terzo: l'autovettore va rinormalizzato a somma 1, perché NumPy lo restituisce di norma euclidea 1 — un vettore legittimo, ma non una distribuzione di probabilità. Il terzo metodo parte dall'uniforme e moltiplica mille volte. I tre coincidono a 1e-10, e la stazionaria è esattamente (6, 6, 4, 4, 2, 3)/25: numeri razionali, perché le probabilità di partenza lo erano.

Che cosa ci dice. La pagina / e la pagina /blog si dividono quasi la metà del traffico, e /strumenti ne prende l'8% pur essendo linkata dalla home esattamente come le altre due: la differenza è che da /strumenti si esce subito e non ci si torna, mentre nella coppia /blog/a1-/blog/a2 si rimbalza. L'ultima riga è il teorema di Kac: il tempo medio di ritorno in uno stato è 1/πᵢ, quindi il crawler ripassa dalla home in media ogni 4,17 clic. È anche il modo giusto di leggere la stazionaria: non «la probabilità finale», ma la frazione di tempo che il processo passa su ciascuna pagina, su un orizzonte lungo.

Il gap spettrale

Sapere che la catena converge non dice quando. La risposta sta nel secondo autovalore. Ordiniamo gli autovalori per modulo decrescente, 1 = |λ₁| > |λ₂| ≥ |λ₃| ≥ …, e supponiamo per semplicità che P sia diagonalizzabile. Allora ogni distribuzione iniziale si scrive come π₀ = π + Σ_{k≥2} cₖ wₖ, dove i wₖ sono gli autovettori sinistri degli autovalori λₖ. Applicando P una volta ciascun pezzo viene moltiplicato per il suo autovalore, e applicandola n volte

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

Il termine π resta fermo perché il suo autovalore è 1; tutti gli altri si contraggono, e quello che muore per ultimo è il più grande in modulo. Quindi ‖πₙ − π‖ ≈ C · |λ₂|ⁿ: l'errore decade geometricamente con ragione |λ₂|. La quantità 1 − |λ₂| si chiama gap spettrale, e il numero di passi necessari perché l'errore si riduca di un fattore e è circa

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

Per scendere sotto una soglia ε servono invece circa ln(1/ε) / (1 − |λ₂|) passi: il gap dà la scala, il logaritmo della precisione dà il fattore moltiplicativo. Vediamo i numeri della nostra catena.

moduli = sorted(np.abs(np.linalg.eigvals(P)), reverse=True)
print("moduli degli autovalori:", "  ".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} passi")

Output ottenuto:

moduli degli autovalori: 1.0000  0.7547  0.6667  0.5000  0.3323  0.3323
|λ₂| = 0.7547 | gap = 1 − |λ₂| = 0.2453 | t_mix ≈ 4.08 passi

Che cosa ci dice. Un solo autovalore ha modulo 1 (coerente con Perron-Frobenius: la catena è irriducibile e aperiodica), gli ultimi due sono una coppia complessa coniugata — perfettamente lecita, perché P non è simmetrica — e il secondo in modulo vale 0,7547. Il gap è ampio, il tempo di mescolamento è di quattro passi: coerente con quello che avevamo visto, cioè che dopo quattro clic la distribuzione aveva già quasi dimenticato il punto di partenza.

Misurato: due siti con la stessa stazionaria

Il gap non si legge dalla stazionaria: due catene possono avere la stessa distribuzione limite e metterci tempi diversissimi ad arrivarci. Costruiamo due mini-siti di sei pagine, entrambi perfettamente simmetrici e quindi entrambi con stazionaria uniforme a 1/6. Nel primo le pagine sono divise in due sezioni da tre, e il navigatore resta nella sua sezione il 99% delle volte. Nel secondo ogni pagina linka tutte le altre, come un menu globale.

def sito_a_sezioni(interno):
    """Sei pagine in due sezioni da tre: con probabilità 'interno' si resta nella sezione."""
    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          # le due pagine sorelle
            else:
                M[i, j] = (1.0 - interno) / 3.0  # le tre pagine dell'altra sezione
    return M

lento = sito_a_sezioni(0.99)
veloce = np.full((6, 6), 1 / 5) - np.eye(6) / 5  # ogni pagina linka tutte le altre

for nome, M in (("sezioni chiuse", lento), ("menu globale", veloce)):
    moduli = sorted(np.abs(np.linalg.eigvals(M)), reverse=True)
    pi = stazionaria_sistema_lineare(M)
    print(f"{nome:>15} | π₁ = {pi[0]:.4f}, uniforme: {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   sezioni chiuse   menu globale")
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 (("sezioni chiuse", lento), ("menu globale", 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} | passi per scendere sotto 0.01: {passi}")

Output ottenuto:

 sezioni chiuse | π₁ = 0.1667, uniforme: True | |λ₂| = 0.9800 | t_mix ≈ 50.00
   menu globale | π₁ = 0.1667, uniforme: True | |λ₂| = 0.2000 | t_mix ≈ 1.25
   n   sezioni chiuse   menu globale
   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
 sezioni chiuse | passi per scendere sotto 0.01: 228
   menu globale | passi per scendere sotto 0.01: 4

Cosa fa il codice. sito_a_sezioni(0.99) costruisce la catena «due sezioni quasi separate»: da ogni pagina, il 99% della probabilità si divide fra le due sorelle di sezione e l'1% fra le tre pagine dell'altra. veloce è il grafo completo con pesi uniformi. Entrambe hanno la stessa stazionaria uniforme — lo verifica np.allclose(pi, 1/6) — ma i secondi autovalori sono 0,98 e 0,20. La tabella misura la distanza L¹ fra la distribuzione al passo n e l'uniforme, partendo da una pagina sola: per il sito a sezioni chiuse la distanza decade esattamente come 0,98ⁿ, e dopo 100 passi è ancora a 0,13.

Che cosa ci dice. Il rapporto fra i due tempi di mescolamento misurati è 228/4 = 57, mentre il rapporto fra le stime 1/(1 − |λ₂|) è 50/1,25 = 40: stesso ordine di grandezza, non lo stesso numero. Lo scarto non è rumore, è il fattore ln(1/ε) che la stima grezza ignora: con ε = 0,01 vale 4,6, e 4,6 × 50 = 230, che è praticamente il 228 misurato. Trade-off. Il gap spettrale è la stima più economica del tempo di convergenza — un solo autovalore — ma è una stima asintotica: descrive il comportamento quando n è grande, e può essere molto ottimista nei primi passi se la costante C davanti è grande. Quando il numero di stati è grande, calcolare λ₂ costa più del simulare: si usa scipy.sparse.linalg.eigs con k=2 e non numpy.linalg.eig, che richiederebbe la matrice densa.

Bilancio dettagliato e reversibilità

C'è una condizione più forte della stazionarietà, e molto più facile da verificare. Una distribuzione π soddisfa il bilancio dettagliato rispetto a P se, per ogni coppia di stati,

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

In parole: il flusso di probabilità che va da i a j pareggia esattamente quello che torna indietro, coppia per coppia. La stazionarietà chiede molto meno: che pareggino i totali in entrata e in uscita di ogni stato, non i singoli scambi. Una catena che soddisfa il bilancio dettagliato si dice reversibile: guardata all'indietro ha la stessa legge di quando la si guarda in avanti.

Il bilancio dettagliato implica la stazionarietà. Somma la relazione su tutti i possibili i: Σᵢ πᵢPᵢⱼ = Σᵢ πⱼPⱼᵢ = πⱼ Σᵢ Pⱼᵢ = πⱼ · 1 = πⱼ. Il membro di sinistra è esattamente la componente j-esima di πP, quindi πP = π. ∎ Il viceversa è falso, e il nostro grafo di crawl lo dimostra: ha una stazionaria, ma il flusso / → /strumenti vale 0,24 × 1/3 = 0,08 mentre il flusso di ritorno /strumenti → / vale 0,08 × 1/2 = 0,04. Pareggiano i conti solo a giro completo, passando per /contatti.

def scarto_bilancio_dettagliato(P, pi):
    """max |πᵢPᵢⱼ − πⱼPⱼᵢ| su tutte le coppie: zero se e solo se la catena è reversibile."""
    F = pi[:, None] * P                     # F[i, j] = flusso di probabilità da i a j
    return np.abs(F - F.T).max()

pi = stazionaria_sistema_lineare(P)
scarto = scarto_bilancio_dettagliato(P, pi)
print(f"grafo di crawl: scarto massimo = {scarto:.6f} | reversibile: {bool(scarto < 1e-12)}")

# stesso sito con i link reciproci: passeggiata su un grafo NON orientato
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("gradi:", [int(g) for g in gradi], "| π = grado/2m :", "  ".join(f"{v:.6f}" for v in pi_r))
print("π R = π :", bool(np.allclose(pi_r @ R, pi_r)))
print("bilancio dettagliato verificato a 1e-12:", bool(scarto_bilancio_dettagliato(R, pi_r) < 1e-12))

Output ottenuto:

grafo di crawl: scarto massimo = 0.040000 | reversibile: False
gradi: [3, 3, 2, 2, 2, 2] | π = grado/2m : 0.214286  0.214286  0.142857  0.142857  0.142857  0.142857
π R = π : True
bilancio dettagliato verificato a 1e-12: True

Che cosa ci dice. Se i link fossero tutti reciproci — cioè se il grafo fosse non orientato — la catena diventerebbe reversibile e la sua stazionaria si leggerebbe senza risolvere nessun sistema: πᵢ = dᵢ / Σⱼ dⱼ, il grado del nodo diviso la somma dei gradi. La verifica è immediata: πᵢ Pᵢⱼ = (dᵢ/2m)(1/dᵢ) = 1/2m, che è simmetrico in i e j. Una formula chiusa al posto di un solve: è questo che si compra con la reversibilità.

Perché conta davvero. Il bilancio dettagliato si usa al contrario. Nella lezione sul campionamento e i metodi Monte Carlo il problema era estrarre campioni da una distribuzione che si sa calcolare solo a meno di una costante. La soluzione è costruire una catena di Markov che abbia quella distribuzione come stazionaria, farla girare, e prendere gli stati visitati come campioni. Progettare una P che soddisfi π = πP direttamente sarebbe difficile; imporre il bilancio dettagliato, che è una condizione locale coppia per coppia, è facile — ed è esattamente ciò che fa il rapporto di accettazione di Metropolis-Hastings. La catena si può progettare per avere la stazionaria che si vuole, e questa equazione è lo strumento con cui la si progetta.

L'esperimento: un navigatore vero, il pozzo e il teletrasporto

Verifica per simulazione: un milione di clic

Tutto quello che abbiamo calcolato finora è algebra lineare: π è la soluzione di un sistema, e nessun navigatore si è mai mosso. La stazionaria però ha un significato operativo preciso — è la frazione di tempo che un singolo navigatore, lasciato camminare abbastanza a lungo, passa su ciascuna pagina. È una promessa verificabile, e la verifichiamo: un navigatore, un milione di clic, seme fissato.

from bisect import bisect_right

def simula(P, partenza, passi, seme=42):
    """Segue un singolo navigatore e conta quante volte visita ogni pagina."""
    rng = np.random.default_rng(seme)
    cumulate = [list(riga) for riga in P.cumsum(axis=1)]   # somme cumulate per riga
    for riga in cumulate:
        riga[-1] = 1.0                                     # 1/3+1/3+1/3 fa 0.9999999999999999
    sorteggi = rng.random(passi)                           # tutti i numeri casuali in un colpo
    conteggi = [0] * len(P)
    stato = partenza
    for u in sorteggi:
        stato = bisect_right(cumulate[stato], u)           # salta da sé gli stati a probabilità 0
        conteggi[stato] += 1
    return np.array(conteggi) / passi

frequenze = simula(P, pagine.index("/contatti"), 1_000_000)
pi = stazionaria_sistema_lineare(P)
print("     pagina   frequenza   stazionaria   scarto")
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"scarto massimo su 10⁶ passi: {np.abs(frequenze - pi).max():.5f}")

Output ottenuto (seme 42, NumPy 2.3.5):

     pagina   frequenza   stazionaria   scarto
          /     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
scarto massimo su 10⁶ passi: 0.00043

Cosa fa il codice. Il trucco per non impiegare mezz'ora è tutto nelle prime righe: le somme cumulate di ogni riga si calcolano una volta sola, i numeri casuali si estraggono tutti insieme, e la scelta dello stato successivo diventa una ricerca binaria su sei numeri. bisect_right fa anche una cosa gratis: sugli stati a probabilità zero la cumulata non cresce, quindi non possono mai essere scelti, e non serve nessun caso speciale. Il seme è fissato a 42, quindi la simulazione è riproducibile cifra per cifra.

Che cosa ci dice. La frequenza di visita misurata e la stazionaria calcolata coincidono con uno scarto massimo di 0.00043. È l'ordine di grandezza giusto: la deviazione standard della frequenza empirica su N passi scala come 1/√N, e 1/√(10⁶) = 0,001. Per guadagnare una cifra decimale servirebbero cento milioni di passi — la stessa legge dei rendimenti decrescenti che governa ogni stima Monte Carlo, e che hai visto nella lezione sul campionamento. Trade-off. Su sei stati la simulazione è un lusso: il sistema lineare dà la risposta esatta in microsecondi. Su milioni di stati si inverte tutto, perché il solve diventa impraticabile e la simulazione resta lineare nel numero di passi: è per questo che nel modulo 08 i ritorni si stimano camminando, non risolvendo.

Il problema del pozzo

Finora il mini-sito era sano. Rompiamolo nel modo più comune in un crawl reale: una pagina senza link uscenti. Succede continuamente — un PDF, una pagina di conferma, una scheda prodotto esaurita. Il crawler ci arriva e non ha dove andare. Modelliamolo come un auto-anello, che è anche quello che fa un browser lasciato aperto: lo stato resta quello.

link_pozzo = dict(link)
link_pozzo["/contatti"] = ["/contatti"]       # nessun link uscente: il crawler resta lì
P_pozzo = matrice_di_transizione(pagine, link_pozzo)

descrivi("con il pozzo", P_pozzo)
pi_pozzo = stazionaria_potenze(P_pozzo, 2000)
print("stazionaria:", "  ".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} passi")

Output ottenuto:

con il pozzo: 2 classe/i, irriducibile = False
   transitoria | ['/', '/blog', '/blog/a1', '/blog/a2', '/strumenti']
    ricorrente | ['/contatti']
stazionaria: 0.0000  0.0000  0.0000  0.0000  0.0000  1.0000
|λ₂| = 0.9278 | t_mix ≈ 13.84 passi

Che cosa ci dice. La catena si è spezzata in due classi: le cinque pagine normali sono diventate transitorie (il crawler ci passa un numero finito di volte, poi le abbandona per sempre) e /contatti è l'unica classe chiusa. La stazionaria esiste ancora ed è unica, ma è degenere: tutta la massa su una pagina, zero su tutte le altre. Come misura di importanza non vale niente — e si noti che il gap spettrale non segnala il problema, anzi: con |λ₂| = 0,9278 la catena converge lentamente, dà numeri che sembrano in movimento, e arriva comunque a una risposta inutile. Un gap piccolo non è un errore; una classe chiusa degenere sì.

La correzione classica è il fattore di teletrasporto di Page, Brin, Motwani e Winograd (1999): con probabilità d il navigatore segue un link, con probabilità 1 − d si annoia e salta su una pagina a caso del sito. La nuova matrice è

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

dove J è la matrice di tutti uno e n il numero di pagine. Tutte le caselle di G sono strettamente positive, quindi la catena è irriducibile e aperiodica per costruzione, e Perron-Frobenius garantisce una stazionaria unica e positiva. In più c'è un risultato esatto sugli autovalori (Langville e Meyer, 2004): se P ha autovalori 1, λ₂, λ₃, …, allora G ha autovalori 1, dλ₂, dλ₃, …. Il teletrasporto non sposta l'autovalore 1 e moltiplica tutti gli altri per d: il gap spettrale si allarga sempre.

def con_teletrasporto(P, d=0.85):
    """G = d·P + (1 − d)/n · J : con probabilità 1 − d si salta su una pagina a caso."""
    n = len(P)
    return d * P + (1 - d) / n * np.ones((n, n))

G = con_teletrasporto(P_pozzo)
descrivi("con teletrasporto", 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"|λ₂| prima = {prima:.4f} | dopo = {dopo:.4f} | rapporto = {dopo / prima:.4f}")

Output ottenuto:

con teletrasporto: 1 classe/i, irriducibile = True
    ricorrente | ['/', '/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
|λ₂| prima = 0.9278 | dopo = 0.7886 | rapporto = 0.8500

Cosa fa il codice. con_teletrasporto mescola la matrice originale con la matrice uniforme, con peso d = 0,85 — il valore usato nell'articolo originale di PageRank. La catena torna irriducibile, la stazionaria torna informativa, e il rapporto fra i secondi autovalori è 0,8500, cioè esattamente d: il risultato di Langville e Meyer, verificato numericamente fino alla quarta cifra. Il tempo di mescolamento scende da 13,84 a 4,73 passi, e questo spiega perché PageRank si calcoli con il metodo delle potenze: con d = 0,85 bastano poche decine di iterazioni, su qualunque grafo.

Il limite della correzione, dichiarato. Guarda il numero di /contatti: 0,4765, quasi metà di tutta la massa. Il teletrasporto ha riparato la matematica — esistenza, unicità, convergenza — ma non ha tolto al pozzo il suo vantaggio strutturale: ci si entra da due parti e non se ne esce mai se non per noia. Il rimedio completo, nell'articolo del 1999, è trattare le pagine senza link uscenti ridistribuendo la loro massa su tutte le pagine, invece di lasciarle ciclare su sé stesse. La regola da portare via. Il teletrasporto è una correzione del modello, non una pezza numerica: cambia la catena di cui stai calcolando la stazionaria, e chi legge il risultato ha il diritto di sapere che d = 0,85 è una scelta tua.

Dal processo alla decisione

Finora il processo si subiva. Il crawler non sceglie: la matrice P è data, e lui la segue. Scrivi la transizione come P(s′ | s) — la probabilità di finire nello stato s′ partendo da s — e hai descritto tutto quello che può succedere.

Aggiungi adesso un'unica cosa: in ogni stato puoi compiere un'azione, scelta in un insieme finito, e l'azione cambia la matrice di transizione. La transizione diventa

P(s′ | s, a),

cioè non più una matrice ma una matrice per ogni azione. Questo oggetto si chiama processo decisionale di Markov, e una regola che a ogni stato associa un'azione si chiama politica: fissata una politica, le matrici tornano a essere una sola, e sei di nuovo in una catena di Markov, con la sua stazionaria e il suo gap spettrale. È per questo che la catena viene prima: un processo decisionale è una famiglia di catene, indicizzata dalle politiche, più un criterio per sceglierne una.

Di qui in avanti il modulo 08 aggiunge due ingredienti che qui non servono: una ricompensa associata alle transizioni e un criterio per confrontare politiche su un orizzonte lungo. Da lì nasce l'equazione di Bellman. Ma l'oggetto su cui tutto questo poggia — lo spazio degli stati, la transizione probabilistica, la proprietà che il futuro dipende solo dal presente — è quello che hai costruito in questa lezione, e il vocabolario è lo stesso: stato, transizione, distribuzione, convergenza.

Errori tipici (sintomo, causa, rimedio)

1. Convenzione riga/colonna scambiata. Sintomo: π non torna, oppure «converge» a una distribuzione che non c'entra niente con il grafo. Causa: hai scritto P @ pi con P stocastica per righe, cioè hai simulato la catena trasposta. Rimedio: dichiara la convenzione in un commento, verifica P.sum(axis=1) e non P.sum(axis=0), e controlla il primo passo a mano su uno stato con una sola uscita.
2. Cercare una convergenza che non esiste. Sintomo: il metodo delle potenze oscilla fra due o più valori e non si ferma mai. Causa: la catena è periodica, e Perron-Frobenius garantisce la convergenza solo con l'aperiodicità. Rimedio: calcola il periodo, o guarda quanti autovalori hanno modulo 1; se ne hai più di uno, la stazionaria esiste ma non è un limite, e va calcolata risolvendo il sistema lineare.
3. Autovettore non normalizzato a somma 1. Sintomo: le componenti di π sommano a 1,7 o a −1. Causa: numpy.linalg.eig restituisce autovettori di norma euclidea 1, con il segno deciso dall'algoritmo. Rimedio: dividi per la somma, non per la norma, e se la somma è negativa cambia segno a tutto il vettore prima di dividere.
4. Autovettore destro al posto del sinistro. Sintomo: ottieni un vettore di componenti tutte uguali, che sembra una distribuzione uniforme plausibile. Causa: hai chiamato eig(P) invece di eig(P.T), e l'autovettore destro di autovalore 1 è sempre il vettore di soli uni. Rimedio: se il risultato è uniforme su una catena che non è simmetrica, hai preso il lato sbagliato.
5. Trattare come markoviano uno stato che non lo è. Sintomo: il modello prevede male, e peggiora dove la storia conta (code, ritardi, penalità). Causa: lo stato è troppo povero. Rimedio: elenca le informazioni che useresti per fare la previsione a mano; se ne usi una che non è nello stato, lo stato va allargato — e dichiara quanto cresce lo spazio.
6. Metodo delle potenze fermato troppo presto. Sintomo: due esecuzioni con numero di iterazioni diverso danno risposte diverse alla terza cifra. Causa: il gap spettrale è piccolo, e dopo n passi l'errore è ancora |λ₂|ⁿ. Rimedio: non fissare le iterazioni, fissa una soglia sullo scarto fra due passi consecutivi, e come controllo calcola |λ₂| una volta sola per sapere che ordine di grandezza aspettarti.
7. Leggere la stazionaria come «la risposta finale». Sintomo: si conclude che il sistema «finirà» in quello stato. Causa: confusione fra distribuzione limite e destinazione. Rimedio: la stazionaria è la frazione di tempo passata in ciascuno stato sul lungo periodo — il processo continua a muoversi, sempre.

I sette errori si riconoscono dal sintomo: quattro danno numeri sbagliati ma plausibili (convenzione, autovettore destro, normalizzazione, potenze fermate presto) e tre sono errori di modello (periodicità ignorata, stato troppo povero, stazionaria letta come destinazione). Contro i primi quattro basta calcolare π in due modi e confrontare, contro i tre di modello serve guardare il grafo prima delle formule.

Best practice

• Dichiara la convenzione riga o colonna nella prima riga di commento, e verificala con un assert sulle somme: costa una riga e previene l'errore più costoso.
• Controlla irriducibilità e periodo prima di cercare una stazionaria, non dopo: sono due visite del grafo e un gcd, e dicono se la domanda che stai facendo ha una risposta.
• Calcola π in due modi indipendenti e confronta: sistema lineare e metodo delle potenze usano algoritmi diversi, e un disaccordo è quasi sempre un errore di costruzione della matrice.
• Guarda il secondo autovalore prima di scegliere il numero di iterazioni, e usa una soglia sullo scarto invece di un numero fisso di passi.
• Dichiara i parametri del modello — il fattore di teletrasporto, l'ipotesi di uniformità sui link, l'omogeneità nel tempo — accanto al risultato: sono scelte tue, non proprietà dei dati.
• Per matrici grandi e sparse usa scipy.sparse.linalg.eigs con k piccolo, e non costruire mai Pⁿ: moltiplica un vettore per volta.

Le sei regole sono controlli di costo molto diverso: la convenzione e i due modi di calcolo costano una riga ciascuno, il controllo di irriducibilità e periodo due visite del grafo, la dichiarazione dei parametri è una frase. Se devi sceglierne una sola, scegli il confronto fra due metodi indipendenti: smaschera quasi tutti gli errori di costruzione della matrice.

Glossario minimo

• Proprietà di Markov (Markov property): il futuro dipende dal passato solo attraverso lo stato presente.
• Matrice di transizione (transition matrix): Pᵢⱼ = probabilità di passare da i a j; righe non negative che sommano a 1.
• Irriducibile (irreducible): da ogni stato si raggiunge ogni altro stato; una sola classe di comunicazione.
• Aperiodica (aperiodic): il massimo comun divisore dei tempi di ritorno vale 1; senza questa proprietà la distribuzione oscilla.
• Ricorrente / transitorio (recurrent / transient): si torna con probabilità 1 / si torna solo un numero finito di volte.
• Distribuzione stazionaria (stationary distribution): π = πP; autovettore sinistro di autovalore 1, normalizzato a somma 1.
• Gap spettrale (spectral gap): 1 − |λ₂|; governa la velocità con cui l'errore |λ₂|ⁿ si annulla.
• Tempo di mescolamento (mixing time): numero di passi per avvicinarsi alla stazionaria entro una soglia; circa ln(1/ε)/(1 − |λ₂|).
• Bilancio dettagliato (detailed balance): πᵢPᵢⱼ = πⱼPⱼᵢ; implica la stazionarietà, non viceversa.
• Teorema di Kac: il tempo medio di ritorno in uno stato è 1/πᵢ.

Le dieci voci seguono l'ordine delle sezioni e si dividono in tre gruppi: la catena (proprietà, matrice), la struttura del grafo (irriducibile, aperiodica, ricorrente e transitorio) e la dinamica (stazionaria, gap spettrale, mescolamento, bilancio dettagliato, Kac). Dopo la definizione di stazionaria le altre voci sono sue conseguenze o suoi strumenti di misura.

Mini-progetto

Costruisci un modulo che, dato un dizionario di link, risponda a tutte le domande della lezione.

1. Scrivi analizza(link) che restituisce: la matrice di transizione, le classi di comunicazione con l'etichetta ricorrente/transitoria, il periodo, la stazionaria calcolata in due modi indipendenti e |λ₂|;
2. applicalo a tre grafi: quello della lezione, lo stesso con /contatti senza link uscenti, e lo stesso con il teletrasporto a d = 0,85; stampa per ciascuno le cinque risposte in una tabella;
3. gestisci e documenta i due casi limite: una pagina che linka sé stessa e basta (classe chiusa di un solo stato) e un grafo con due componenti separate (stazionaria non unica: il tuo codice deve accorgersene e dirlo, invece di restituire un numero a caso).

Le tre richieste costruiscono un unico modulo che deve rispondere sempre alle stesse cinque domande su qualunque grafo: la prima definisce l'interfaccia, la seconda lo mette alla prova su tre casi, la terza lo spinge ai casi limite dove le ipotesi di Perron-Frobenius cadono. I casi limite sono la parte che distingue uno script da un modulo.

Consegna: uno script che stampa la tabella delle tre catene e, in fondo, una riga che dichiara quale ipotesi del modello consideri meno difendibile e perché.

Challenge

• Easy. Carta e penna: data P = [[0.9, 0.1], [0.5, 0.5]], calcola π risolvendo π = πP con π₁ + π₂ = 1. Soluzione: la prima equazione è π₁ = 0,9π₁ + 0,5π₂, cioè 0,1π₁ = 0,5π₂, cioè π₁ = 5π₂. Con il vincolo 5π₂ + π₂ = 1 si ottiene π₂ = 1/6 e π₁ = 5/6. Verifica con la seconda equazione: 0,1·(5/6) + 0,5·(1/6) = 1/12 + 1/12 = 1/6 ✓.
• Medium. Carta e penna: dimostra che il bilancio dettagliato implica la stazionarietà, e poi costruisci un controesempio al viceversa con tre stati (suggerimento: un ciclo orientato 1 → 2 → 3 → 1 con probabilità tutte uguali a 1 ha stazionaria uniforme, ma i flussi in un verso non pareggiano quelli nell'altro).
• Hard. Codice: costruisci due catene a quattro stati con la stessa stazionaria uniforme, una periodica di periodo 2 e una aperiodica, e mostra numericamente la differenza: per la prima il metodo delle potenze oscilla, per la seconda converge; conta gli autovalori di modulo 1 in entrambi i casi.
• Expert. Codice: stima il tempo di mescolamento per simulazione — parti da ogni stato possibile, misura dopo quanti passi la distanza L¹ dalla stazionaria scende sotto 0,01, prendi il massimo sugli stati di partenza — e confrontalo con 1/(1 − |λ₂|) e con ln(100)/(1 − |λ₂|) su almeno cinque catene diverse. Quale delle due stime è sistematicamente migliore?
• Enterprise. Prendi uno stato di un sistema del tuo lavoro che evolve a passi discreti — gli stati di un ticket, di un ordine, di un lead nel funnel — e stima P dai dati storici contando le transizioni (P̂ᵢⱼ = nᵢⱼ / Σₖ nᵢₖ). Riporta la stazionaria, il tempo medio di permanenza in ciascuno stato e |λ₂|. Poi la parte difficile: verifica se la proprietà di Markov regge davvero, confrontando la distribuzione delle transizioni uscenti da uno stato in funzione dello stato precedente. Se le due distribuzioni sono diverse, lo stato va allargato, e il tuo modello a un passo sta mentendo.

I cinque livelli salgono dal conto a mano a un'analisi su dati tuoi: i primi due si risolvono con carta e penna (la stazionaria di una matrice 2×2 e il legame fra bilancio dettagliato e stazionarietà), il terzo e il quarto chiedono codice per confrontare catene periodiche e aperiodiche e stime del tempo di mescolamento, l'ultimo stima P dai tuoi dati e verifica che la proprietà di Markov regga davvero.

Gli esercizi di codice di questa lezione sono in un solo file, esercizi.py (circa 11 KB). Quattro esercizi di codice sui temi della lezione: costruire la matrice di transizione da un dizionario di link, calcolare la stazionaria risolvendo il sistema lineare, verificarla simulando un navigatore per un milione di passi, e stimare il tempo di mescolamento confrontandolo con il gap spettrale. Scheletro con le parti da completare e, in fondo, le soluzioni commentate. Richiede Python 3.12+ e NumPy 2.3. Eseguilo con python3 esercizi.py: gira anche con le parti da completare lasciate vuote e ti dice quali mancano.

Recap

Mappa mentale. Proprietà di Markov: il futuro dipende dal passato solo attraverso lo stato — e lo stato lo scegli tu. P stocastica per righe, distribuzione come vettore riga, πₙ = π₀Pⁿ. Irriducibile ⇒ la stazionaria esiste, è unica, è positiva. Aperiodica ⇒ ci si arriva. π = πP ⟺ Pᵀπᵀ = πᵀ: autovettore sinistro di autovalore 1. L'errore decade come |λ₂|ⁿ, il tempo di mescolamento è circa ln(1/ε)/(1 − |λ₂|). Bilancio dettagliato πᵢPᵢⱼ = πⱼPⱼᵢ ⇒ stazionarietà, non viceversa. Una catena si subisce; con P(s′ | s, a) si decide.

Cheat sheet.

3 concetti: lo stato riassume il passato, non lo cancella; la stazionaria è l'autovettore sinistro di autovalore 1; il secondo autovalore governa la velocità. 3 errori: convenzione riga/colonna scambiata; cercare un limite su una catena periodica; usare l'autovettore destro al posto del sinistro. 3 best practice: verificare irriducibilità e periodo prima di calcolare π; calcolarla in due modi indipendenti; fermare il metodo delle potenze su una soglia, non su un numero fisso di passi. 1 esercizio: il mini-progetto di questa lezione.

Prossima lezione: finisce qui la matematica. Dalla prossima lezione i modelli: il processo di un progetto di machine learning, e le tre metriche che mentono quando le classi sono sbilanciate.

Domande per andare avanti

• Come si stima una matrice di transizione dai dati quando alcune transizioni non si sono mai osservate, e come si evita di dichiararle impossibili? (Correzione di Laplace e prior di Dirichlet: lezione 13 di questo modulo.)
• Che cosa succede quando lo spazio degli stati è continuo o infinito, e la matrice diventa un nucleo di transizione? (È il caso dei metodi Monte Carlo a catena di Markov: lezione 14 di questo modulo.)
• Come si sceglie l'azione migliore quando a ogni stato ne hai molte e il risultato si vede solo dopo molti passi? (Processo decisionale di Markov ed equazione di Bellman, modulo 08.)

Le tre domande sono ganci verso altre lezioni e ciascuna ha già una risposta tra parentesi: la prima riguarda le transizioni mai osservate e rimanda alla correzione di Laplace, la seconda gli spazi degli stati continui e i metodi Monte Carlo a catena di Markov, la terza la scelta delle azioni nel processo decisionale di Markov. Servono a decidere che cosa leggere dopo.