Algoritmi randomizzati - Monte Carlo, Las Vegas e amplificazione

Probabilità di base: Spazi di probabilità discreti e uniformiIn uno spazio discreto la probabilità è determinata dalla densità discreta p(ω) = P({ω}), con somma 1, e P(A) è la somma di p(ω) sugli esiti di A; negli spazi uniformi (esiti equiprobabili) P(A) = |A| / |Ω|, casi favorevoli su casi possibili.Spazi di probabilità discreti e uniformi →, Indipendenza di eventiA e B sono indipendenti se P(A ∩ B) = P(A) P(B), cioè se sapere che uno si è verificato non cambia la probabilità dell'altro; l'indipendenza passa ai complementari, non va confusa con l'incompatibilità, e per più eventi va richiesta su ogni sottofamiglia.Indipendenza di eventi →, Valore attesoIl valore atteso E[X] = Σ x p_X(x) è la media dei valori di X pesata con le loro probabilità (esiste se la serie converge assolutamente); per una funzione g vale E[g(X)] = Σ g(x) p_X(x) senza trovare la legge di g(X), ed E è lineare: E[aX + bY + c] = aE[X] + bE[Y] + c.Valore atteso →, Distribuzione geometricaGeo(p) è il numero della prova in cui arriva il primo successo in prove indipendenti: P(X = n) = (1−p)^(n−1) p per n ≥ 1, P(X > n) = (1−p)^n (lunga attesa), media 1/p, varianza (1−p)/p², ed è senza memoria.Distribuzione geometrica →. Complessità: Complessità in tempo e caso pessimoPerché lo studio sperimentale non basta; complessità al caso pessimo come massimo sul numero di operazioni tra le istanze di una data taglia; stima con limiti superiore e inferiore senza trovare l'istanza peggiore; esempi arrayMax, prefixAverages e InsertionSort; efficienza asintotica e limiti dell'analisi.Complessità in tempo e caso pessimo →, Notazione asintoticaDefinizioni di O, Omega, Theta e o piccolo con le costanti c ed n0; esempi con costanti esplicite; proprietà (polinomi, esponenziali, logaritmi, somme, implicazioni tra notazioni); sommatorie notevoli; terminologia (logaritmica, lineare, polinomiale, esponenziale).Notazione asintotica →.

1. Il modello

Un algoritmo randomizzato è un algoritmo che, oltre all'input, legge un flusso di bit casuali indipendenti ed equi (lanci di moneta). Il suo comportamento è quindi una variabile aleatoria, definita sullo spazio dei bit casuali.

Garanzia. Le probabilità sono sui bit casuali, non sull'input. Un enunciato come «per ogni input $x$, $\Pr[\text{errore}]\le 1/4$» vale nel caso peggiore sull'input: nessun avversario che scelga $x$ in anticipo può fare peggio. È diverso dall'analisi del caso medio, dove si assume una distribuzione sugli input. Un esempio piccolo in cui una scelta casuale batte ogni strategia deterministica nel caso pessimo: Esercizio - Due scatole di caramelle e soglia casuale contro ogni input (lezioni MIT).

Perché servono: sono spesso più semplici e più veloci di quelli deterministici (verifica di prodotti di matrici, test di primalità, taglio minimo: Union bound e verifica randomizzata - identità polinomiali e prodotti di matriciSpazi discreti richiamati; inclusione-esclusione (con le troncate di Bonferroni) e union bound Pr[unione] ≤ somma, valido senza indipendenza e stretto quando gli eventi sono rari (collisioni: C(n,2)/m). Principio delle decisioni differite: si fissa tutto tranne una scelta casuale e si limita la probabilità su quella. Lemma di Schwartz-Zippel: un polinomio non nullo di grado totale d, valutato in un punto di S^n scelto uniformemente, si annulla con prob ≤ d/|S| (induzione sul numero di variabili); test di identità polinomiale con errore unilaterale, ripetuto (d/|S|)^k. Freivalds: per verificare AB=C si controlla A(Br)=Cr con r in {0,1}^n in tempo O(n^2), errore ≤ 1/2 (stretto), 2^-k con k prove. Fingerprint: due stringhe di n bit si confrontano con a mod p per un primo casuale p ≤ M, errore ≤ n/π(M).Union bound e verifica randomizzata - identità polinomiali e prodotti di matrici →, Taglio minimo randomizzato - algoritmo di KargerTaglio minimo globale di un multigrafo con n vertici. Algoritmo di Karger: finché restano più di 2 vertici, scegliere un arco uniformemente tra i rimasti e contrarlo (fondere gli estremi, tenere i paralleli, togliere i cappi); gli archi tra i due supervertici finali sono un taglio. Un taglio minimo C di k archi sopravvive se nessun suo arco viene contratto; con n_i vertici il grado minimo è ≥ k, quindi |E_i| ≥ n_i k/2 e Pr[contrarre un arco di C] ≤ 2/n_i; il prodotto telescopico dà Pr[output = C] ≥ 2/(n(n-1)) = 1/C(n,2). Con t = C(n,2) ln(1/δ) ripetizioni indipendenti il fallimento è ≤ δ. Corollario: un grafo ha al più C(n,2) tagli minimi (stretto per il ciclo). Karger-Stein: contrarre fino a n/√2+1 vertici (sopravvive con prob ≥ 1/2) e ricorrere due volte, T(n) = O(n² log n), successo Ω(1/log n).Taglio minimo randomizzato - algoritmo di Karger →). Un esempio già noto è quicksort con pivot casuale, che diventa Las Vegas (Algoritmi di ordinamentoSelection sort e insertion sort (quadratici, con invarianti), merge sort (divide et impera, Theta(n log n)), quick sort (caso pessimo quadratico, medio n log n), heap sort (in loco, Theta(n log n)), ordinamento senza confronti per chiavi intere in un intervallo piccolo, limite inferiore Omega(n log n) per gli algoritmi basati su confronti.Algoritmi di ordinamento →, Analisi probabilistica di quicksort e quickselectCon pivot casuale due elementi di rango $i<j$ vengono confrontati con probabilità $2/(j-i+1)$; sommando sulle coppie il numero atteso di confronti di quicksort è esattamente $2(n+1)H_n-4n\approx2n\ln n-2{,}85,n$ (stessa formula dalla ricorrenza $C_n=n-1+\frac2n\sum_{k<n}C_k$). Quickselect trova il $k$-esimo elemento con attesa lineare, al più $4n$ confronti (mediana $\approx3{,}39,n$). Pivot fisso su input casuale equivale a pivot casuale su input qualunque.Analisi probabilistica di quicksort e quickselect →).

2. Las Vegas e Monte Carlo

Definizione. Un algoritmo è Las Vegas se la risposta è sempre corretta e il tempo di esecuzione $T$ è una variabile aleatoria: si studia $E[T]$. È Monte Carlo se il tempo è limitato (deterministico) ma la risposta può essere sbagliata con una probabilità limitata.

Esempio (cerca un 1). Un array $A[0..n-1]$ ha almeno metà delle celle uguali a $1$. Si vuole un indice $i$ con $A[i]=1$.

• Las Vegas: si sceglie un indice a caso, si legge, e si ripete finché si trova un $1$. Ogni tentativo riesce con probabilità $p\ge 1/2$ e i tentativi sono indipendenti, quindi $T$ è geometrica e $E[T]=1/p\le 2$.
• Monte Carlo: si fanno $k$ tentativi e, se nessuno riesce, si restituisce l'ultimo indice estratto (sbagliato). L'errore è $(1-p)^k\le 2^{-k}$.

Esecuzione a mano su $A=[0,1,1,0,0,1,0,1]$ ($n=8$, quattro $1$, $p=1/2$) con la sequenza casuale $3,0,6,5,\dots$: $A[3]=0$, $A[0]=0$, $A[6]=0$, $A[5]=1$. Si restituisce $5$ dopo $T=4$ tentativi. In generale $\Pr[T\ge 4]=(1/2)^3=1/8$ (i primi tre falliscono), $\Pr[T=4]=1/16$.

import random

def cerca_uno_lv(A):                       # Las Vegas: sempre corretto, tempo aleatorio
    t = 0
    while True:
        t += 1
        i = random.randrange(len(A))
        if A[i] == 1:
            return i, t

def cerca_uno_mc(A, k):                    # Monte Carlo: k tentativi, errore <= 2^-k
    for _ in range(k):
        i = random.randrange(len(A))
        if A[i] == 1:
            return i
    return i                               # sbagliato con prob (1-p)^k

A = [0, 1, 1, 0, 0, 1, 0, 1]
N = 100000
print(sum(cerca_uno_lv(A)[1] for _ in range(N)) / N)                 # circa 2
print(sum(A[cerca_uno_mc(A, 3)] == 0 for _ in range(N)) / N)         # circa 1/8

3. Errore unilaterale e bilaterale, classi di complessità

Per un problema di decisione (linguaggio $L$: si risponde SI se $x\in L$, NO altrimenti) si distinguono tre casi di Monte Carlo.

Definizione (errore unilaterale). Un algoritmo è RP (randomized polynomial time) se lavora in tempo polinomiale e

• $x\in L\ \Rightarrow\ \Pr[\text{SI}]\ge 1/2$;
• $x\notin L\ \Rightarrow\ \Pr[\text{SI}]=0$.

Un SI è quindi certo, un NO può essere un errore. La classe co-RP è speculare: un NO è certo, un SI può essere un errore.

Definizione (errore bilaterale). Un algoritmo è BPP (bounded-error probabilistic polynomial time) se, in tempo polinomiale, risponde correttamente con probabilità almeno $2/3$ per ogni $x$, sia in $L$ sia fuori.

Definizione. ZPP (zero-error) è la classe dei linguaggi decisi da un algoritmo Las Vegas con tempo atteso polinomiale.

In breve vale $\text{P}\subseteq\text{ZPP}\subseteq\text{RP}\subseteq\text{BPP}$. Le costanti $1/2$ e $2/3$ sono arbitrarie, perché l'amplificazione (sezione 4) le porta ovunque in $(0,1)$ e $(1/2,1)$. Per esempio $\text{RP}\subseteq\text{BPP}$ perché un errore unilaterale $\le1/2$, ripetuto due volte, diventa $\le1/4<1/3$.

Teorema. $\text{ZPP}=\text{RP}\cap\text{co-RP}$.

Dimostrazione. ($\supseteq$) Siano $A$ un algoritmo RP e $B$ un co-RP per $L$. A turno: si esegue $A$ e, se dice SI, si risponde SI (certo); si esegue $B$ e, se dice NO, si risponde NO (certo); altrimenti si ricomincia. Se $x\in L$, $A$ dice SI con probabilità $\ge 1/2$ (e $B$ non dice mai NO); se $x\notin L$, $B$ dice NO con probabilità $\ge 1/2$. Ogni giro termina quindi con probabilità $\ge 1/2$: i giri sono geometrici, in media $\le 2$, e la risposta è sempre corretta. ($\subseteq$) Sia $V$ un Las Vegas con $E[T]\le q(n)$ polinomiale. Lo si interrompe a $2q(n)$ passi: per Markov $\Pr[T>2q(n)]\le 1/2$. Se si risponde NO allo scadere del tempo si ottiene un algoritmo RP (un SI è sempre certo, perché prodotto da $V$ concluso, e se $x\in L$ si risponde SI con probabilità $\ge1/2$); se si risponde SI si ottiene un co-RP. $\square$

4. Amplificazione

Errore unilaterale: ripetizione

Si esegue l'algoritmo $k$ volte con bit casuali indipendenti e si risponde SI se almeno una esecuzione dice SI.

• Se $x\notin L$ nessuna esecuzione dice SI: errore $0$.
• Se $x\in L$, ogni esecuzione sbaglia con probabilità $\le 1-p$, e per indipendenza tutte sbagliano con probabilità $(1-p)^k$.

Poiché $1-x\le e^{-x}$ (la retta tangente in $0$ sta sotto la convessa $e^{-x}$), $$\Pr[\text{errore}]\le(1-p)^k\le e^{-pk}\le\varepsilon\qquad\text{se}\qquad k\ge\frac{\ln(1/\varepsilon)}{p}.$$

Esempi numerici (verificati): con $p=1/2$ e $\varepsilon=2^{-20}$ servono esattamente $k=20$ ripetizioni (la formula con $e$ dice $28$); con $p=0{,}3$ e $\varepsilon=10^{-6}$ servono $k=\lceil\ln10^{-6}/\ln0{,}7\rceil=39$ (la formula con $e$ dice $47$, pessimistica ma semplice). Il costo cresce linearmente in $k$, l'errore cala esponenzialmente: per questo gli algoritmi randomizzati sono pratici.

Errore bilaterale: voto di maggioranza

Con errore bilaterale non si può fidare di un SI isolato. Si eseguono $k$ esecuzioni indipendenti ($k$ dispari) e si risponde come la maggioranza. Ogni esecuzione è giusta con probabilità $p\ge2/3$, quindi il numero di esecuzioni giuste è $X\sim\text{Bin}(k,p)$ (Bernoulli, binomiale e uniforme discretaBernoulli Be(p): un solo tentativo, vale 1 con probabilità p; binomiale Bin(n, p): numero di successi in n prove indipendenti, P(X = k) = C(n,k) p^k (1−p)^(n−k), media np e varianza np(1−p); uniforme discreta: n valori equiprobabili.Bernoulli, binomiale e uniforme discreta →) e $$\Pr[\text{errore}]=\Pr\Big[X\le\tfrac{k-1}{2}\Big]=\sum_{j=0}^{(k-1)/2}\binom kj p^j(1-p)^{k-j}.$$

Valori esatti per $p=2/3$ (calcolati con Python):

Per $k=3$: $\binom30(1/3)^3+\binom31(2/3)(1/3)^2=\frac1{27}+\frac6{27}=\frac7{27}=0{,}2593$, cioè sbaglia se le giuste sono $0$ o $1$ su $3$.

Quanto decresce? Con Chebyshev (Disuguaglianze di Markov, Chebyshev e JensenMarkov: per X ≥ 0, P(X ≥ a) ≤ E[X]/a; Chebyshev: P(|X − μ| ≥ ε) ≤ Var(X)/ε²; Jensen: per φ convessa, φ(E[X]) ≤ E[φ(X)]. Stimano probabilità e medie conoscendo solo media e varianza.Disuguaglianze di Markov, Chebyshev e Jensen →): $E[X]=2k/3$, $\text{Var}[X]=2k/9$, l'evento $X\le k/2$ richiede uno scarto di almeno $k/6$, quindi $\Pr\le\frac{2k/9}{(k/6)^2}=\frac8k$. È solo $1/k$: per $k=81$ dà $0{,}099$, contro il valore vero $0{,}00097$. Il limite di Chernoff-Hoeffding (Limiti di Chernoff e di HoeffdingMarkov applicato a e^(tX) dà code esponenziali. Per somme di Bernoulli indipendenti con media μ: P(X ≥ (1+δ)μ) ≤ [e^δ/(1+δ)^(1+δ)]^μ ≤ e^(-δ²μ/3) per 0<δ≤1, P(X ≤ (1-δ)μ) ≤ e^(-δ²μ/2), P(X ≥ R) ≤ 2^(-R) se R ≥ 6μ. Hoeffding per X_i in [a_i,b_i]: P(|S-E[S]| ≥ t) ≤ 2exp(-2t²/Σ(b_i-a_i)²). Applicazioni: n ≥ ln(2/δ)/(2ε²) campioni per stimare una probabilità, amplificazione BPP con k = O(log 1/ε) ripetizioni, bilanciamento di insiemi con scarto √(4m ln n).Limiti di Chernoff e di Hoeffding →) dà invece $$\Pr[\text{errore}]\le e^{-2k(p-1/2)^2}=e^{-k/18}\qquad(p=2/3),$$ esponenziale in $k$. Per ottenere errore $\le0{,}01$ servono $47$ esecuzioni esatte (la formula di Hoeffding ne chiede $\lceil18\ln100\rceil=83$). Serve però uno scarto costante da $1/2$: se $p=1/2+1/\text{poly}(n)$ il numero di ripetizioni sale a $\text{poly}(n)$.

5. Da Las Vegas a Monte Carlo e ritorno

Las Vegas $\to$ Monte Carlo (Markov). Sia $\mu=E[T]$. Per la disuguaglianza di Markov $$\Pr[T\ge c\mu]\le\frac1c.$$ Si interrompe l'algoritmo dopo $c\mu$ passi e, se non ha finito, si restituisce una risposta qualsiasi: si ottiene un Monte Carlo con errore $\le 1/c$ e tempo $c\mu$. Ripetendo $k$ blocchi indipendenti di durata $c\mu$ e tenendo la prima risposta completata, l'errore cala a $c^{-k}$.

Nell'esempio «cerca un 1» $\mu=2$: fermarsi a $10$ tentativi ($c=5$) dà per Markov errore $\le0{,}2$, mentre il valore vero è $\Pr[T>10]=2^{-10}\approx0{,}00098$. Markov sfrutta solo la media, quindi è largo ma sempre valido.

Monte Carlo $\to$ Las Vegas. Si può fare solo se la risposta si può verificare in fretta: si ripete il Monte Carlo finché il verificatore accetta. Con errore $\varepsilon$ il numero di ripetizioni è geometrico, in media $1/(1-\varepsilon)$. Nell'esempio, il verificatore legge $A[i]$. Senza un verificatore (per esempio un Monte Carlo bilaterale su un problema di decisione) non si sa in generale farlo: che BPP sia uguale a ZPP è una questione aperta.

6. Esempio: test di primalità di Miller-Rabin

Dato $n$ dispari si scrive $n-1=2^sd$ con $d$ dispari. Un numero $a\in[1,n-1]$ è un testimone di composizione se $$a^d\not\equiv1\pmod n\quad\text{e}\quad a^{2^rd}\not\equiv-1\pmod n\ \text{ per }r=0,1,\dots,s-1.$$

• Se $n$ è primo non esistono testimoni: in $\mathbb Z_n$ (un campo, CampiUn campo è un insieme di numeri con somma e prodotto che rispettano le regole usuali (associativa, commutativa, neutri, opposti, inversi, distributiva): è ciò che serve per risolvere le equazioni di primo grado. Esempi: Q, R, C e i campi finiti.Campi →) l'unica radice quadrata di $1$ oltre a $1$ è $-1$, e la catena $a^d,a^{2d},\dots,a^{n-1}=1$ (piccolo teorema di Fermat) deve partire da $1$ oppure passare per $-1$.
• Se $n$ è composto almeno i $3/4$ dei valori $a$ sono testimoni (teorema di Rabin, senza dimostrazione qui).

L'algoritmo sceglie $a$ a caso e risponde COMPOSTO se è un testimone, PROBABILMENTE PRIMO altrimenti. L'errore è unilaterale (un COMPOSTO è certo) e $\le1/4$ per esecuzione: con $k$ esecuzioni indipendenti $\le4^{-k}$, per la formula della sezione 4 con $p=3/4$. Il valore $4^{-k}$ è la probabilità che un composto fissato superi il test; la probabilità che un numero estratto a caso e promosso dal test sia composto si ricava con la regola di Bayes e dipende da quanti primi ci sono (Esercizio - Moneta a due teste e regola di Bayes per le risposte di un test randomizzato (lezioni MIT)).

A mano, $n=221=13\cdot17$. $n-1=220=2^2\cdot55$, quindi $s=2$, $d=55$.

• $a=2$: $2^{55}\bmod221=128$ (non è $1$ né $220$) e $128^2=16384\equiv30$ (non è $220$). Nessuna condizione di non-testimone è soddisfatta: $2$ è un testimone, quindi $221$ è composto.
• $a=174$: $174^{55}\equiv47$ e $47^2=2209\equiv220=-1$: qui $a$ è un «bugiardo» (passa il test pur con $n$ composto).

Contando tutte le basi con Python: per $n=221$ i bugiardi sono solo $\{1,21,47,174,200,220\}$, cioè $6$ su $220$ ($97{,}3\%$ di testimoni, ben oltre $3/4$). Il test di Fermat da solo ($a^{n-1}\equiv1$) è molto peggiore: su $561=3\cdot11\cdot17$ (primo numero di Carmichael) tutte le $320$ basi coprime con $561$ lo superano, mentre Miller-Rabin ha solo $10$ bugiardi su $560$.

import random

def testimone(a, n):                       # True se a prova che n e' composto
    d, s = n - 1, 0
    while d % 2 == 0:
        d //= 2; s += 1
    x = pow(a, d, n)
    if x in (1, n - 1):
        return False
    for _ in range(s - 1):
        x = x * x % n
        if x == n - 1:
            return False
    return True

def miller_rabin(n, k=20):                 # n dispari >= 5; errore <= 4^-k
    return not any(testimone(random.randrange(2, n - 1), n) for _ in range(k))

print(miller_rabin(221), miller_rabin(561), miller_rabin(1000003))   # False False True

Domande d'esame

1. Differenza tra Las Vegas e Monte Carlo, con un esempio. Las Vegas: risposta sempre giusta, tempo aleatorio (quicksort con pivot casuale, «cerca un 1»). Monte Carlo: tempo fisso, risposta sbagliata con probabilità limitata (Miller-Rabin).
2. Dimostra che $\text{ZPP}=\text{RP}\cap\text{co-RP}$. Vedi il teorema della sezione 3: i due algoritmi a turno nel verso $\supseteq$, Markov a $2q(n)$ passi nel verso $\subseteq$.
3. Un algoritmo RP ha $p=1/2$. Quante ripetizioni per errore $\le2^{-30}$? $(1/2)^k\le2^{-30}$ dà $k=30$; la formula con $e$ ne chiede $\lceil30\ln2/0{,}5\rceil=42$.
4. Un Las Vegas ha tempo atteso $100$ ms. Come ottenere un Monte Carlo con errore $\le1\%$? Markov con $c=100$: interromperlo a $10$ s; l'errore è $\le1/100$.