Abbiamo visto nel precedente articolo che un algoritmo Monte Carlo serve a calcolare una media della forma seguente:
$$\langle O \rangle = \sum_{\{s\}} O(\{s\}) p(\{s\})$$
dove gli $\{s\}$ sono una configurazione di spin, $O$ è un osservabile e $p$ è una distribuzione di probabilità difficilmente calcolabile analiticamente, nella pratica $p = \frac{1}{Z}e^{-\beta H}$. L’idea è di produrre una serie di configurazioni di spin
$$\{s^{(1)}\}, \{s^{(2)}\}, \ldots, \{s^{(T)}\} \tag{1}$$
secondo la probabilità $p(\{s\})$, dove cioè le configurazioni appaiono con frequenza relativa data dalla loro probabilità. Una volta campionata in questa maniera la distribuzione di probabilità, possiamo calcolare la media di un’osservabile $O$ tramite
$$\langle O \rangle \approx \frac{1}{T}\sum_{i=1}^T O(\{s^{(i)}\})$$
Monte Carlo permette di produrre i campioni $(1)$ senza calcolare esplicitamente la probabilità $p({s})$, ed è quindi vantaggioso se il calcolo della probabilità è computazionalmente costoso. Qui vediamo esplicitamente il funzionamento di un algoritmo particolare detto algoritmo di Metropolis.
L’algoritmo di Metropolis
Vediamo come produrre la serie di configurazioni $(1)$ secondo l’algoritmo Metropolis. Partiamo da una configurazione $\{s^{(1)}\}$ qualsiasi: possiamo scegliere tutti gli spin uguali, o tutti a caso, oppure in qualsiasi altra maniera. Il risultato finale non dipende dalla scelta iniziale di configurazione. Ad ogni passaggio applichiamo una certa operazione $K$ e otteniamo la configurazione successiva:
$$\{s^{(i+1)}\} = K( \{s^{(i)}\} )$$
La serie così ottenuta è esattamente il campionamento $(1)$ che cerchiamo se $K$ è scelta in maniera appropriata. Se vogliamo aumentare il numero $T$ di configurazioni generate, basta aumentare il numero di passaggi dell’algoritmo, che può andare avanti anche all’infinito.
L’operazione $K$ che svolgiamo ad ogni passaggio dipende solo dalla configurazione attuale. Al momento $i$, la configurazione è $\{s^{(i)}\}$ ed effettuiamo le seguenti operazioni:
- Consideriamo tutti gli spin nel reticolo. Per ogni spin, effettuiamo quanto segue:
- Al tempo $i$, nella configurazione attuale $\{s^{(i)}\}$, lo spin selezionato ha valore $s$. Scegliamo a caso un valore $s’$ dello spin diverso da $s$ e chiamiamo $\{s^{\prime (i)}\}$ la configurazione in cui abbiamo cambiato $s \to s’$ per lo spin in questione, mentre tutti gli altri sono immutati.
- Calcoliamo il rapporto tra le probabilità dopo il cambio e prima del cambio, $$r = \frac{p(\{s^{\prime (i)}\})}{p(\{s^{(i)}\})}$$
- Se $r \geq 1$, allora effettuiamo il cambio $s\to s’$ per lo spin in questione; altrimenti effettuiamo il cambio con probabilità $r$.
- Torniamo al punto $2$ finché non abbiamo esaurito tutti gli spin sul reticolo. Alla fine del processo abbiamo ottenuto la configurazione $\{s^{(i+1)}\}$.
Il vantaggio è tutto nel punto $3$. In un sistema di fisica statistica, la probabilità è data da
$$p(\{s\}) = \frac{1}{Z}e^{-\beta H(\{s\})}$$
dove $H$ è l’Hamiltoniana. Pertanto il rapporto tra le due probabilità è dato da
$$r = \exp{\bqty{-\beta \pqty{ H(\{s^{\prime (i)}\}) -H(\{s^{(i)}\}) }}}$$
Quindi se $r \geq 1$, cioè l’energia diminuisce, accettiamo sempre il cambio; altrimenti lo accettiamo con probabilità $r$. In particolare, notiamo che in $r$ la funzione di partizione è stata semplificata, e quindi non c’è alcun bisogno di calcolarla. Inoltre tipicamente l’Hamiltoniana è locale, cioè ogni spin è collegato solo a quelli vicini. Quindi poiché cambiamo un solo spin per volta, nella differenza tra le due Hamiltoniane la maggior parte dei termini è immutata e quindi si cancella. Il calcolo di $r$ è perciò estremamente più semplice del calcolo di $p$: qui sta tutta l’efficienza di Monte Carlo.
Tipicamente, le prime configurazioni generate non hanno ancora “termalizzato“, cioè non hanno la corretta distribuzione, e quindi vanno scartate. Vedremo questo fenomeno più nel dettaglio nei prossimi articoli.
Nel prossimo articolo vedremo perché questo metodo funziona e poi nei successivi un esempio pratico col modello di Ising.