Catene di Markov: la matrice che si applica per sempre#
Per mescolare bene un mazzo di cinquantadue carte servono sette mescolate. Lo dimostra un risultato di matematica pura, al quale nel gennaio del 1990 il New York Times dedicò un articolo. Dave Bayer e Persi Diaconis avevano studiato che cosa succede a un mazzo mescolato alla maniera dei croupier, dividendolo in due e facendo cadere le carte a intreccio, e avevano scoperto un fenomeno inatteso.
Fino a quattro mescolate la distanza fra la distribuzione dei mazzi ottenuti e quella di un mazzo in ordine casuale, misurata in variazione totale, resta praticamente uguale a uno. In altri termini, chi conosce la procedura riconosce quasi sempre un mazzo così mescolato. Alla settima la distanza scende per la prima volta sotto la metà, e da lì in poi si dimezza a ogni mescolata [BD92]. Fra il mazzo ancora riconoscibile e quello indistinguibile da uno casuale, quindi, non c’è una salita graduale, c’è una soglia, e sta fra la quinta mescolata e la settima.
Quel risultato riguarda un mazzo di carte, ma la sua struttura è quella di una catena di Markov. Un mazzo mescolato, infatti, è uno stato; la mescolata è una regola che porta da uno stato al successivo, con un po” di caso; e la domanda «dopo quante volte non si distingue più da un mazzo qualsiasi?» riguarda la velocità con cui la catena si avvicina all’equilibrio. Per quasi tutte le catene il secondo autovalore della matrice che descrive la regola ne governa la velocità; il mazzo è uno dei casi in cui quel numero da solo non basta.
Il patto: domani dipende solo da oggi#
Nel gioco dell’oca, dove finisci al prossimo turno dipende da due cose: la casella su cui ti trovi adesso e il numero che esce dal dado. Non conta come ci sei arrivato, se hai fatto un percorso fortunato o disastroso, se hai già perso tre turni. Tutta la tua storia è riassunta da un solo dato, la casella.
Questo è il patto di Markov, e i sistemi che lo rispettano si chiamano catene di Markov. È una restrizione seria, e proprio per questo porta lontano. Un sistema che ricorda tutto è impossibile da studiare, mentre uno che ricorda solo dove si trova si maneggia con l’algebra lineare.
Un accorgimento torna utile in continuazione. Se il gioco avesse una regola del tipo «chi è già finito in prigione una volta, la seconda volta esce subito», il patto sarebbe rotto: la casella da sola non basterebbe più a decidere il futuro. Basta però cambiare che cosa si chiama stato. Invece della sola casella si prende la coppia (casella, sei già stato in prigione?), e il patto torna valido, con il doppio degli stati. In generale, quasi ogni sistema si può rendere markoviano allargando la definizione di stato, ma gli stati si moltiplicano.
Descritto il gioco, si vuole sapere dove si finisce. Non in che casella cadrai tu stasera, che dipende dai dadi, ma come si distribuiscono i giocatori sulle caselle dopo molti turni: quale frazione sta qui, quale lì. Il gioco parte da una distribuzione (all’inizio tutti sulla casella di partenza) e a ogni turno la rimescola; la domanda è dove quella distribuzione va a finire.
Nei conti la distribuzione è un vettore: una lista di numeri, uno per casella (la frazione di giocatori che sta lì). Le regole del gioco stanno in una tabella di numeri, la matrice di transizione, con una colonna per ogni casella di partenza e dentro le probabilità di arrivare in ciascuna casella. Un turno di gioco è la tabella moltiplicata per la lista. In concreto, per ogni casella d’arrivo si fa uno scontrino, probabilità per quota di giocatori, e si somma.
Un processo \(\{S_t\}_{t\ge 0}\) a valori in un insieme finito \(\{1,\dots,n\}\) è una catena di Markov se
cioè il futuro è condizionatamente indipendente dal passato, dato il presente. I numeri \(P_{ij}\) formano la matrice di transizione \(\mathbf{P}\in\mathbb{R}^{n\times n}\).
I testi si dividono su una convenzione, e noi prendiamo \(\mathbf{P}\) stocastica per colonne: \(P_{ij}\) è la probabilità di andare in \(i\) partendo da \(j\), quindi ogni colonna è una distribuzione di probabilità e somma a \(1\). Con questa scelta lo stato del sistema è un vettore colonna \(\mathbf{x}_t\in\mathbb{R}^n\) di probabilità, e l’evoluzione è la moltiplicazione a sinistra:
Molti testi di probabilità adottano la convenzione trasposta (righe stocastiche, vettori riga, \(\boldsymbol{\pi}_{t+1} = \boldsymbol{\pi}_t \mathbf{P}\)), e confrontando una formula con la sua fonte conviene guardare da che parte sta la somma a uno.
L’ipotesi di Markov è meno restrittiva di quanto sembri, per via dell’allargamento dello stato. Un processo che dipende dagli ultimi \(k\) istanti, infatti, diventa markoviano prendendo come stato la \(k\)-upla degli ultimi \(k\) valori, al costo di \(n^k\) stati. È la costruzione dei modelli linguistici che condizionano ogni parola alle \(k\) precedenti (gli \(n\)-gram, dove \(n=k+1\) conta le parole della finestra e non gli stati), e la stessa con cui un problema di decisione in cui conta la storia recente si riporta a un processo decisionale di Markov.
Dove si finisce: un autovettore di autovalore uno#
Fig. 3.25 La distribuzione \(\mathbf{x}_t\) su quattro pagine del web (è la catena del PageRank, costruita nel paragrafo «Il web come catena: PageRank»), partendo da quattro valori uguali: le colonne sono le probabilità di trovarsi su ciascuna pagina, le righe tratteggiate i valori della distribuzione stazionaria. Le colonne A, B e C oscillano intorno al proprio valore con scarti sempre più piccoli; la D, che nessuna pagina linka, lo raggiunge già dopo un passo e non cambia più.#
Applicare la regola una volta cambia la distribuzione. Applicarla di nuovo la cambia ancora, ma di meno. Andando avanti si arriva quasi sempre a una distribuzione che la regola lascia com’è: applicandola, non succede più niente. Quella si chiama distribuzione stazionaria, e Fig. 3.25 la mostra come la riga tratteggiata verso cui le colonne si assestano.
Il disegno mostra la fotografia d’insieme, che a un certo punto smette di cambiare. I singoli giocatori, invece, continuano a cambiare casella a ogni turno: la fotografia resta ferma perché quanti entrano in una casella e quanti ne escono si pareggiano.
La distribuzione stazionaria, però, si può trovare risolvendo, invece che aspettando. Cercare una distribuzione che la regola lascia identica a sé stessa equivale a cercare una direzione che la matrice non devia, cioè un autovettore con fattore di allungamento uguale a uno. È una delle venature del legno di una matrice; qui la venatura che conta è quella che non si allunga né si accorcia, e il sistema ci scivola sopra e ci resta.
Una distribuzione \(\boldsymbol{\pi}\in\mathbb{R}^n\), con \(\pi_i\ge 0\) e \(\sum_i\pi_i=1\), si dice stazionaria se
il che la rende un autovettore di \(\mathbf{P}\) relativo all’autovalore \(1\). L’autovalore \(1\) esiste sempre. Le colonne di \(\mathbf{P}\), infatti, sommano a \(1\), quindi \(\mathbf{1}^\top\mathbf{P} = \mathbf{1}^\top\), e \(\mathbf{1}\) è autovettore sinistro con autovalore \(1\); poiché \(\mathbf{P}\) e \(\mathbf{P}^\top\) hanno lo stesso spettro, esiste anche l’autovettore destro.
Osserviamo però che questo ragionamento non garantisce che l’autovettore si possa scegliere a componenti non negative, e quindi normalizzare a una distribuzione: lo garantisce la parte del teorema di Perron-Frobenius che vale per ogni matrice a entrate non negative, e senza di essa l’autovettore potrebbe avere segni misti.
Vale inoltre \(|\lambda|\le 1\) per ogni autovalore. Gli autovettori non sono distribuzioni, quindi l’argomento deve valere per vettori qualsiasi: una matrice stocastica non allunga la norma \(\ell_1\) di nessun vettore, dato che \(\sum_i \big|\sum_j P_{ij} x_j\big| \le \sum_j |x_j| \sum_i P_{ij} = \sum_j |x_j|\), e un autovettore con \(|\lambda|>1\) la violerebbe. Lo spettro sta quindi nel disco unitario, con almeno un punto sul bordo.
La stazionaria si calcola in tre modi, tutti usati:
decomposizione spettrale, prendendo l’autovettore di autovalore \(1\) e normalizzandolo a somma uno;
sistema lineare \((\mathbf{P}-\mathbf{I})\boldsymbol{\pi}=\mathbf{0}\) con il vincolo \(\mathbf{1}^\top\boldsymbol{\pi}=1\), un sistema indeterminato reso determinato dalla normalizzazione, nella forma della sezione sui sistemi lineari;
iterazione della potenza, applicando \(\mathbf{P}\) ripetutamente a una distribuzione qualsiasi. È il metodo che si usa sulle matrici enormi, dove fattorizzare è fuori discussione e moltiplicare per una matrice sparsa costa poco.
L’equilibrio è dinamico: la catena continua a saltare da uno stato all’altro, mentre il flusso complessivo si pareggia. Una condizione più forte, il bilancio dettagliato \(P_{ij}\pi_j = P_{ji}\pi_i\), chiede che si pareggi ogni singola coppia di stati. È sufficiente per la stazionarietà, e lo si vede sommando su \(j\): \(\sum_j P_{ij}\pi_j=\pi_i\sum_j P_{ji}=\pi_i\). Non è necessaria, ma è comoda, perché per costruire una catena con una stazionaria assegnata basta imporre una condizione locale fra coppie di stati, senza risolvere nessun sistema. È la proprietà su cui si costruiscono i campionatori usati nelle macchine di Boltzmann.
Una catena che soddisfa il bilancio dettagliato si dice reversibile, e ne segue un fatto spettrale: posto \(\mathbf{D}=\operatorname{diag}(\boldsymbol{\pi})\), con \(\pi_i>0\), la matrice \(\mathbf{D}^{-1/2}\mathbf{P}\mathbf{D}^{1/2}\) è simmetrica, quindi gli autovalori di \(\mathbf{P}\) sono tutti reali. Per questo gli autovalori complessi della catena del PageRank, e la rotazione che producono, sono la firma di una catena non reversibile.
Quando la distribuzione stazionaria è una sola#
Per una catena finita una distribuzione stazionaria esiste sempre. È unica se la catena è irriducibile (da ogni stato si può arrivare a ogni altro). Se la catena è anche aperiodica (non ha un ritmo fisso), la distribuzione al tempo \(t\) converge a essa da qualunque partenza. Le due condizioni si capiscono meglio guardando come si rompono.
Prima condizione: da ogni stato si deve poter arrivare a ogni altro. Se il tabellone è fatto di due zone senza passaggi fra loro, chi parte nella prima ci resta per sempre, e le distribuzioni di equilibrio sono almeno due: una per zona. Chiedere «dove finisce il sistema» non ha una risposta sola, perché dipende da dove è partito. È il caso di un web fatto di due gruppi di pagine che non si citano a vicenda, una situazione più comune di quanto si creda.
Seconda condizione: non deve esserci un ritmo fisso. Due stanze sono collegate da una porta girevole che a ogni passo obbliga a cambiare stanza. Partendo dalla prima sarai nella seconda a ogni turno dispari e nella prima a ogni turno pari, per sempre: la distribuzione oscilla fra due valori e non si assesta mai, anche se la media sul lungo periodo è metà e metà.
Basta una sola possibilità di restare fermi, o due giri di ritorno le cui lunghezze non abbiano divisori in comune (uno di due passi e uno di tre, per esempio), per rompere il ritmo e far convergere tutto. Un solo giro di tre passi, invece, non basta: la distribuzione gira fra tre valori, come fra due nella porta girevole.
Con tutte e due le condizioni la garanzia è forte: esiste una sola distribuzione di equilibrio, e la catena ci arriva da qualunque punto di partenza. La memoria dell’inizio si consuma, e l’equilibrio dipende soltanto dalla regola. È il motivo per cui un mazzo mescolato abbastanza a lungo non ricorda più com’era ordinato all’inizio.
Una catena è irriducibile se per ogni coppia \((i,j)\) esiste \(t\) con \((\mathbf{P}^t)_{ij}>0\) (il grafo delle transizioni è fortemente connesso), ed è aperiodica se il massimo comun divisore dei tempi \(t\ge 1\) con \((\mathbf{P}^t)_{ii}>0\) (quelli in cui un ritorno è possibile) vale \(1\).
Teorema di Perron-Frobenius, nella forma che serve qui: se \(\mathbf{P}\) è stocastica, irriducibile e aperiodica, allora l’autovalore \(1\) è semplice, tutti gli altri hanno modulo strettamente minore di \(1\), l’autovettore associato si può scegliere a componenti tutte positive e risulta unico a meno della normalizzazione. Di conseguenza
I due controesempi mostrano che nessuna delle due ipotesi è di comodo. Con \(\mathbf{P}=\mathbf{I}\) (riducibile all’estremo) ogni distribuzione è stazionaria e non c’è unicità. Con \(\mathbf{P}=\begin{pmatrix}0&1\\1&0\end{pmatrix}\) (irriducibile ma periodica di periodo \(2\)) lo spettro è \(\{1,-1\}\): la stazionaria \((\tfrac12,\tfrac12)\) esiste ed è unica, ma \(\mathbf{P}^t\mathbf{x}_0\) oscilla e non converge, perché l’autovalore \(-1\) sta anch’esso sul bordo del disco.
Un’ipotesi più debole basta per l’esistenza e l’unicità senza la convergenza. L’irriducibilità da sola garantisce una stazionaria unica e strettamente positiva, e per di più il teorema ergodico: qualunque sia la partenza, \(\frac1T\sum_{t<T} f(S_t)\to\sum_i\pi_i f(i)\) quasi certamente per ogni funzione \(f\) degli stati, e in particolare la frazione di tempo passata in \(i\) tende a \(\pi_i\). L’aperiodicità serve soltanto alla convergenza della distribuzione al tempo \(t\), non a quella delle medie nel tempo: la catena che salta fra due stati ha le medie giuste anche se la sua distribuzione oscilla. Il teorema ergodico, non la convergenza di \(\mathbf{P}^t\mathbf{x}_0\), giustifica un campionatore MCMC, che stima un valore atteso con la media lungo una sola traiettoria.
La catena pigra \(\tfrac12(\mathbf{I}+\mathbf{P})\) si usa soprattutto nelle dimostrazioni: sposta lo spettro nel disco di centro \(\tfrac12\) e raggio \(\tfrac12\), che tocca il bordo solo in \(1\), al costo di dimezzare circa il gap.
Quanto ci si mette: il secondo autovalore#
Sapere che la catena converge non basta, perché la regola si applica un numero finito di volte e serve sapere dopo quante. Il modulo del secondo autovalore, \(|\lambda_2|\), dà la velocità asintotica di avvicinamento all’equilibrio; il tempo di mescolamento misura invece quanti passi servono.
Lo scarto fra la distribuzione di adesso e quella di equilibrio, a lungo andare, si accorcia a ogni passo della stessa frazione, almeno in media. È come un libretto di risparmio al contrario, che ogni giorno perde la stessa quota. Ad esempio, se ogni passo dimezza lo scarto, dopo dieci passi ne resta circa un millesimo, perché un mezzo moltiplicato per sé stesso dieci volte fa uno su 1024. Fra due passi vicini il rapporto può anche oscillare, e allora a essere costante è la media su molti passi: è quella la velocità della catena.
Quella frazione è un numero che si legge nella matrice, ed è il secondo dei suoi fattori di allungamento, il più grande dopo l’uno che tiene ferma la distribuzione di equilibrio. Se è piccolo si arriva in fretta, se è vicino a uno ci si mette moltissimo. Poiché lo scarto si moltiplica per quella frazione, il numero di passi non cresce in proporzione alla precisione voluta. In altre parole, ogni cifra decimale in più aggiunge sempre lo stesso numero di passi, quale che sia la cifra.
Il mazzo di carte è un caso speciale. Il ritmo regolare vale comunque a lungo andare, ma prima di arrivarci il passaggio dal mescolato male al mescolato bene non è graduale: per quattro o cinque mescolate lo scarto resta grande, poi crolla in un paio di passaggi, poi non c’è quasi più niente da guadagnare. Certe famiglie di catene enormi si comportano così, con una soglia netta, e questo ha un vantaggio pratico: esiste un numero giusto di passi, e farne il doppio è spreco.
Siano \(1 = |\lambda_1| > |\lambda_2| \ge \dots \ge |\lambda_n|\) gli autovalori di \(\mathbf{P}\). Scomponendo \(\mathbf{x}_0 = \boldsymbol{\pi} + \mathbf{r}_0\) con \(\mathbf{r}_0\) nel sottospazio \(\{\mathbf{r}:\mathbf{1}^\top\mathbf{r}=0\}\), che \(\mathbf{P}\) lascia invariato perché \(\mathbf{1}^\top\mathbf{P}=\mathbf{1}^\top\),
e lo scarto decade geometricamente con ragione \(|\lambda_2|\). La scrittura \(O(|\lambda_2|^t)\) presuppone però che \(\mathbf{P}\) sia diagonalizzabile. Se \(\lambda_2\) ha un blocco di Jordan di taglia \(m\), lo scarto va come \(t^{m-1}|\lambda_2|^t\), con la stessa ragione asintotica ma con un transitorio polinomiale davanti. La quantità
governa il tempo di mescolamento. Con la distanza in variazione totale \(\lVert\boldsymbol{\mu}-\boldsymbol{\nu}\rVert_{TV}=\tfrac12\sum_i|\mu_i-\nu_i|\) e il caso peggiore sulla partenza, si pone
Per una catena reversibile (quella in bilancio dettagliato), irriducibile e aperiodica valgono le due stime [LPW17]
con \(|\lambda_2|\) il più grande modulo fra gli autovalori diversi da \(1\) e \(\pi_{\min}=\min_i\pi_i\). Il fattore \(\log(1/\pi_{\min})\) non è innocuo. Ad esempio, con la stazionaria uniforme sui \(52!\) ordini di un mazzo vale circa \(156\) (il mescolamento a intreccio non è reversibile, e la stima non gli si applica alla lettera, ma l’effetto del numero degli stati è quello). Per una catena non reversibile, come la matrice di Google, lo spettro da solo non basta nemmeno a questa maggiorazione, perché conta anche quanto sono vicini fra loro gli autovettori.
Il logaritmo dice che guadagnare una cifra decimale richiede sempre lo stesso numero di passi, e il fattore davanti dice quanti. Se \(\lambda_2\) è complesso, il decadimento è geometrico soltanto in modulo, con una rotazione sovrapposta. Di conseguenza il rapporto fra due passi consecutivi oscilla, e la stima di \(|\lambda_2|\) si prende come media geometrica su molti passi. Un gap piccolo (una catena quasi riducibile, con due gruppi di stati collegati da pochi passaggi) rende il mescolamento lentissimo, la difficoltà tipica dei campionatori usati in statistica bayesiana.
Il fenomeno del mazzo di carte ha un nome, cutoff (Aldous e Diaconis, 1986): per certe famiglie di catene la distanza dalla stazionaria resta vicina al massimo per un tempo, poi cala a zero in una finestra molto più stretta del tempo stesso. Per il mescolamento a intreccio di \(n\) carte, con \(k=\tfrac32\log_2 n+\theta\) mescolate, la distanza in variazione totale tende a \(1-2\Phi\!\big(-2^{-\theta}/(4\sqrt3)\big)\), con \(\Phi\) la funzione di ripartizione normale [BD92].
Per un mazzo di cinquantadue carte il centro della finestra è \(\tfrac32\log_2 52\approx8{,}55\), dove la distanza vale ancora circa \(0{,}11\); scende sotto \(1/2\) per \(\theta\approx-2{,}2\), quindi a \(k\approx6{,}3\); per questo la settima mescolata è la prima sotto la metà (valori esatti: \(0{,}614\) alla sesta, \(0{,}334\) alla settima). Il decadimento geometrico vale sempre asintoticamente, ma nella finestra del cutoff descrive male ciò che si osserva, e quindi la sola conoscenza di \(|\lambda_2|\) a volte non basta.
Il web come catena: PageRank#
L’esempio che ha reso queste idee note fuori dalla matematica è il modo in cui si ordinano i risultati di una ricerca.
Immagina qualcuno che naviga a caso: apre una pagina, sceglie a caso uno dei link che ci trova, ci va, e ricomincia. Dopo moltissimi salti, la frazione di tempo che passa su ciascuna pagina è la sua importanza. In altre parole, una pagina è importante se molte pagine importanti la linkano. La definizione sembra circolare, ma è quella di una distribuzione stazionaria, e ha una soluzione sola.
Il navigatore così com’è si incastra però in due modi, e le due condizioni dicono come. Può finire su una pagina senza link uscenti e restarci; e può finire in un gruppo di pagine che si linkano solo fra loro e non uscirne più, accumulando tutta l’importanza. La toppa è la stessa per entrambi ed è semplicissima: ogni tanto, con una probabilità piccola, il navigatore si annoia e salta su una pagina qualsiasi del web, scelta a caso (è il teletrasporto). Nella formulazione originale succede circa una volta ogni sette clic.
Questo salto fa tre cose insieme. Rende il web attraversabile da qualunque punto (prima condizione), toglie ogni ritmo fisso (seconda condizione), e accelera la convergenza: più spesso il navigatore si annoia, prima la distribuzione si stabilizza. Il salto, nato per riparare due difetti del modello, è anche ciò che rende possibile il calcolo su miliardi di pagine.
Sia \(\mathbf{M}\) la matrice del web resa stocastica per colonne (\(M_{ij}=1/\ell_j\) se la pagina \(j\) linka la \(i\), con \(\ell_j\) il numero di link uscenti da \(j\)). La matrice di Google è
con \(d\approx 0{,}85\) il fattore di smorzamento e \(n\) il numero di pagine. Il secondo termine è il teletrasporto uniforme; le pagine senza link uscenti si trattano a parte, sostituendo la loro colonna nulla con la distribuzione uniforme.
\(\mathbf{G}\) ha tutte le entrate strettamente positive, quindi è irriducibile e aperiodica, e Perron-Frobenius garantisce che il PageRank \(\mathbf{G}\boldsymbol{\pi}=\boldsymbol{\pi}\) esista, sia unico e abbia componenti positive. Il risultato che rende il calcolo possibile su scala reale è però una stima sul secondo autovalore:
indipendentemente da come è fatto il grafo dei link. Lo smorzamento, introdotto come correzione del modello, si rivela così una garanzia di convergenza a velocità nota.
Il limite si dimostra con lo stesso argomento usato per \(|\lambda|\le1\). Sia \(\lambda\ne1\) un autovalore di \(\mathbf{G}\) con autovettore \(\mathbf{v}\). Da \(\mathbf{1}^\top\mathbf{G}=\mathbf{1}^\top\) segue \(\lambda\,\mathbf{1}^\top\mathbf{v}=\mathbf{1}^\top\mathbf{G}\mathbf{v}=\mathbf{1}^\top\mathbf{v}\), quindi \(\mathbf{1}^\top\mathbf{v}=0\) e il termine di teletrasporto si annulla: \(\mathbf{G}\mathbf{v}=d\,\mathbf{M}\mathbf{v}\). Allora \(|\lambda|\,\lVert\mathbf{v}\rVert_1=d\,\lVert\mathbf{M}\mathbf{v}\rVert_1\le d\,\lVert\mathbf{v}\rVert_1\), perché \(\mathbf{M}\) è stocastica per colonne.
Lo stesso passaggio dà, per ogni scarto \(\mathbf{r}\) con \(\mathbf{1}^\top\mathbf{r}=0\), \(\lVert\mathbf{G}\mathbf{r}\rVert_1\le d\,\lVert\mathbf{r}\rVert_1\). Con \(d=0{,}85\), quindi, ogni iterazione taglia almeno il \(15\%\) dello scarto misurato come somma degli scostamenti (in altre norme il singolo passo può anche peggiorare). Servono circa 43 iterazioni per ridurre lo scarto di un fattore \(10^3\) e 114 per un fattore \(10^8\); ciascuna è un prodotto matrice-vettore sparso, di costo proporzionale al numero di link.
Vale \(|\lambda_2|=d\) se \(\mathbf{M}\) ha almeno due sottoinsiemi chiusi irriducibili [HK03]; nel web di quattro pagine dell’esempio, invece, \(|\lambda_2|=0{,}601<d\).
import numpy as np
# quattro pagine, e chi punta a chi
link = {0: [1, 2], 1: [2], 2: [0], 3: [0, 2]}
nomi, n, d = ["A", "B", "C", "D"], 4, 0.85
M = np.zeros((n, n))
for da, verso in link.items():
for a in verso:
M[a, da] = 1 / len(verso) # colonne stocastiche
print(M.sum(axis=0)) # ogni colonna somma a uno -> [1. 1. 1. 1.]
G = d * M + (1 - d) / n * np.ones((n, n))
# prima strada: l'autovettore di autovalore uno
autoval, autovet = np.linalg.eig(G)
i = np.argmin(np.abs(autoval - 1))
pi = np.real(autovet[:, i]); pi = pi / pi.sum()
print(np.round(pi, 4)) # -> [0.3797 0.1989 0.3839 0.0375]
# seconda strada: applicare la regola finché non si muove più
x = np.ones(n) / n
for passo in range(1, 41):
x = G @ x
print(np.round(x, 4), np.abs(x - pi).max() < 1e-9) # stesso risultato -> True
I due metodi danno lo stesso vettore. Il primo scompone la matrice nei suoi
autovalori e autovettori (eig), il secondo la applica quaranta volte a una
distribuzione uniforme, e dopo quaranta passi lo scarto massimo è sotto
\(10^{-9}\). La pagina \(D\), che nessuna pagina linka, prende il \(3{,}75\%\), pari
a \((1-d)/n\), perché quel poco le arriva soltanto dal teletrasporto.
Il secondo autovalore si legge nella stessa decomposizione, e dice quanto in fretta la seconda strada converge.
print(np.round(autoval, 3))
# -> [ 1. +0.j -0.425+0.425j -0.425-0.425j -0. +0.j ]
secondo = sorted(np.abs(autoval))[-2]
print(round(secondo, 4), d) # -> 0.601 0.85, e vale sempre <= d
# lo scarto si riduce davvero di quel fattore a ogni passo?
x, scarti = np.ones(n) / n, []
for passo in range(40):
x = G @ x
scarti.append(np.abs(x - pi).max())
print([float(round(scarti[k + 1] / scarti[k], 3))
for k in range(14, 19)])
# -> [0.425, 0.5, 0.85, 0.723, 0.425]
print(round((scarti[34] / scarti[14]) ** (1 / 20), 4)) # -> 0.601
Il rapporto fra due scarti consecutivi oscilla con periodo quattro, fra \(0{,}425\) e \(0{,}85\). Lo spiegano i due autovalori \(-0{,}425\pm0{,}425\mathrm{j}\), che sono numeri complessi (in Python \(\mathrm{j}\) è l’unità immaginaria, la radice di \(-1\)).
Moltiplicare per un numero complesso contrae di un fattore, il modulo, e ruota di un angolo, l’argomento. Qui il modulo è \(0{,}601\) e l’angolo tre ottavi di giro, perciò a ogni passo lo scarto si accorcia e gira. Dopo quattro passi ha girato di un giro e mezzo e punta nel verso opposto, e una misura della lunghezza non distingue un vettore dal suo opposto: per questo il ritmo è di quattro passi.
Su venti passi, cioè cinque periodi, la rotazione non cambia più la lunghezza e resta la sola contrazione, tanto che la media geometrica dei rapporti dà \(0{,}601=|\lambda_2|\). In definitiva, la velocità di convergenza si misura su questo inviluppo, mai su un passo solo.
In pratica, con NumPy#
import numpy as np
P = np.array([[0.9, 0.5], # colonne stocastiche: P[i, j] = da j vai in i
[0.1, 0.5]])
# dove si finisce, applicando la regola 50 volte: le due colonne coincidono,
# cioè la partenza non conta più
print(np.round(np.linalg.matrix_power(P, 50), 4))
# -> [[0.8333 0.8333]
# [0.1667 0.1667]]
# fra gli autovalori ce n'è sempre uno che vale 1, ma in virgola mobile può
# uscire 0.9999999999999998: si cerca il più vicino a uno, non l'uguale
w, v = np.linalg.eig(P)
print(np.round(w, 4), np.argmin(np.abs(w - 1))) # -> [1. 0.4] 0
# la stazionaria come sistema lineare: (P - I) pi = 0 con somma uno
A = np.vstack([P - np.eye(2), np.ones(2)])
b = np.array([0.0, 0.0, 1.0])
print(np.round(np.linalg.lstsq(A, b, rcond=None)[0], 4)) # -> [0.8333 0.1667]
Il sistema è sovradeterminato (tre equazioni, due incognite) e si risolve con i minimi quadrati della sezione sulle proiezioni. La riga di normalizzazione aggiunta in fondo sceglie, fra gli infiniti multipli dell’autovettore, l’unico che sia una distribuzione.
Le catene di Markov sotto altri nomi#
I modelli n-gram sono catene di Markov sulle parole, con lo stato allargato alle ultime \(n-1\). I processi decisionali di Markov del reinforcement learning sono catene in cui a ogni passo qualcuno sceglie un’azione; fissata la regola di scelta, il valore di ogni stato risolve un sistema lineare come quelli della sezione sui sistemi lineari[1].
I campionatori delle macchine di Boltzmann costruiscono una catena apposta perché la sua stazionaria sia la distribuzione che si vuole campionare, e ne aspettano il mescolamento. Anche il processo di andata dei modelli di diffusione, che aggiunge rumore un passo alla volta, è una catena di Markov su uno spazio continuo, con una regola che cambia a ogni passo: fuori, quindi, dai teoremi per le catene finite. La sua distribuzione tende comunque a un rumore gaussiano standard, e quei modelli imparano a percorrere la catena al contrario.
Da ricordare
Una catena di Markov è un sistema in cui il prossimo passo dipende solo da dove ti trovi adesso, e non da come ci sei arrivato. Quando la regola sembra ricordare qualcosa di più, quasi sempre basta allargare la definizione di stato e il patto torna valido.
Applicando la regola molte volte la distribuzione si assesta su una configurazione che la regola lascia identica: la distribuzione stazionaria. Il sistema continua a muoversi, è la fotografia d’insieme a fermarsi.
Trovarla è cercare una direzione che la matrice non devia, un autovettore con fattore di allungamento uguale a uno.
La meta è una sola se da ogni stato si può raggiungere ogni altro; e ci si arriva da ovunque solo se in più non c’è un ritmo fisso. Due stanze e una porta girevole che obbliga a cambiare stanza hanno una meta sola, metà e metà, e non ci si assestano mai.
Quanto ci si mette lo dice il secondo fattore di allungamento della matrice: lo scarto dall’equilibrio si riduce, in media, di quella frazione a ogni passo. Nel PageRank quel numero è tenuto sotto controllo dal salto casuale del navigatore, che rende ben definito il calcolo e accelera anche la convergenza.
Da ricordare
\(\Pr[S_{t+1}\mid S_t,\dots,S_0]=\Pr[S_{t+1}\mid S_t]\) definisce la catena; con \(\mathbf{P}\) stocastica per colonne l’evoluzione è \(\mathbf{x}_t=\mathbf{P}^t\mathbf{x}_0\). Un processo con memoria \(k\) diventa markoviano su \(n^k\) stati.
La stazionaria è l’autovettore di autovalore \(1\): \(\mathbf{P}\boldsymbol{\pi}=\boldsymbol{\pi}\). Esiste sempre (le colonne sommano a uno) e tutti gli autovalori stanno nel disco unitario.
Perron-Frobenius: se \(\mathbf{P}\) è irriducibile e aperiodica, l’autovalore \(1\) è semplice, \(\boldsymbol{\pi}\) è unica e positiva, e \(\mathbf{P}^t\mathbf{x}_0\to\boldsymbol{\pi}\) da ogni partenza. \(\mathbf{P}=\mathbf{I}\) rompe la prima ipotesi, \(\begin{pmatrix}0&1\\1&0\end{pmatrix}\) la seconda.
La velocità è il gap spettrale \(1-|\lambda_2|\): lo scarto decade come \(|\lambda_2|^t\) (con un fattore polinomiale se \(\mathbf{P}\) non è diagonalizzabile) e, per una catena reversibile, \(t_{\text{mix}}(\varepsilon)\le\log\big(1/(\varepsilon\,\pi_{\min})\big)/(1-|\lambda_2|)\), dove \(\log(1/\pi_{\min})\) cresce con il numero degli stati; per una catena non reversibile lo spettro da solo non basta. Alcune famiglie mostrano cutoff (un crollo in una finestra stretta), e lì \(|\lambda_2|\) da solo descrive male il transitorio.
PageRank: \(\mathbf{G}=d\mathbf{M}+\frac{1-d}{n}\mathbf{1}\mathbf{1}^\top\) è positiva, quindi la stazionaria esiste ed è unica; e \(|\lambda_2(\mathbf{G})|\le d\) indipendentemente dal grafo, il che rende l’iterazione della potenza praticabile su miliardi di pagine.
Con le catene la probabilità incontra l’algebra lineare: una matrice applicata per sempre a una distribuzione ha un equilibrio che si trova con un autovettore. Restano due strumenti. Il primo misura la sorpresa con un numero solo, e quindi l’errore di un modello: è la teoria dell’informazione. Il secondo dice che cosa accade a questi calcoli quando li esegue una macchina che scrive poche cifre per numero: è l’analisi numerica.