Algoritmi di streaming - campionamento di riserva e conteggio degli elementi distinti

Definizione (modello a flusso). Gli elementi $x_1,x_2,\dots,x_N$ arrivano uno alla volta, si leggono una volta sola e non si possono riguardare. La memoria è molto minore di $N$ (spesso polilogaritmica) e $N$ può essere ignoto. Gli algoritmi sono randomizzati e approssimati: si chiede una stima entro $\varepsilon$ con probabilità almeno $1-\delta$.

È il caso dei contatori di un router, dei log di un server, dei dati dei sensori. Gli strumenti sono quelli già visti: valore atteso e varianza (Linearità del valore atteso e variabili indicatriciLa linearità $E[X+Y]=E[X]+E[Y]$ vale sempre, anche senza indipendenza; con le variabili indicatrici ($E[I_A]=P(A)$) trasforma un conteggio in una somma di probabilità. Applicazioni: attesa del primo successo $1/p$, collezionista di figurine $nH_n$, punti fissi di una permutazione (media $1$), record (media $H_n$), inversioni ($n(n-1)/4$), attesa condizionata, identità di Wald $E[S_N]=E[N]E[X]$.Linearità del valore atteso e variabili indicatrici →, Momenti e deviazioni - Markov, Chebyshev e mediana randomizzataVar(ΣX_i) = ΣVar + 2ΣCov, e basta l'indipendenza a coppie perché la varianza si sommi. Chebyshev dà code polinomiali (collezionista: P(T ≥ nH_n + cn) ≤ π²/(6c²)); il metodo del secondo momento dà P(X=0) ≤ Var(X)/E[X]². Algoritmo della mediana: campione di n^(3/4) elementi, due estremi d e u a ±√n, ordinamento dei soli elementi tra d e u; fallisce con probabilità ≤ n^(-1/4) ed è Las Vegas con 2n + o(n) confronti.Momenti e deviazioni - Markov, Chebyshev e mediana randomizzata →), l'hashing universale (Hashing randomizzato - famiglie universali, filtri di Bloom e fingerprintHashing con catene e costo atteso 1+alfa; famiglie universali (probabilità di collisione al più 1/m) e costruzione h(x)=((ax+b) mod p) mod m con dimostrazione; hashing perfetto a due livelli (spazio lineare, ricerca O(1) nel caso peggiore); filtri di Bloom: falsi positivi (1-e^(-kn/m))^k e k ottimo (m/n) ln 2; fingerprint di stringhe e numero di bit necessari.Hashing randomizzato - famiglie universali, filtri di Bloom e fingerprint →) e l'amplificazione per ripetizione (Algoritmi randomizzati - Monte Carlo, Las Vegas e amplificazioneUn algoritmo randomizzato usa bit casuali indipendenti; la probabilità è sulle scelte dell'algoritmo, mai sull'input (garanzia per ogni input). Las Vegas: sempre corretto, tempo aleatorio (conta il tempo atteso). Monte Carlo: tempo fisso, risposta sbagliata con probabilità limitata; errore unilaterale (RP: accetta solo se x sta nel linguaggio, accetta con prob almeno 1/2) o bilaterale (BPP: giusto con prob almeno 2/3). ZPP = RP ∩ co-RP = Las Vegas con tempo atteso polinomiale. Amplificazione: unilaterale, k ripetizioni indipendenti, errore (1-p)^k ≤ e^(-pk), serve k ≥ ln(1/ε)/p; bilaterale, voto di maggioranza su k esecuzioni, errore binomiale esatto, esponenziale in k (Chernoff, e^(-k/18) per p=2/3) mentre Chebyshev dà solo 8/k. Las Vegas -> Monte Carlo: interrompere a c volte il tempo atteso (Markov, errore ≤ 1/c). Monte Carlo -> Las Vegas: solo se la risposta si verifica in fretta. Esempio: Miller-Rabin, errore ≤ 1/4 per esecuzione.Algoritmi randomizzati - Monte Carlo, Las Vegas e amplificazione →).

Campionamento di riserva

Problema. Estrarre dal flusso un campione di $k$ elementi, uniforme tra i $\binom Nk$ sottoinsiemi possibili, con memoria $k$ e senza conoscere $N$.

Algoritmo (reservoir sampling). Si tengono i primi $k$ elementi. All'arrivo dell'elemento $i>k$: si estrae $j$ uniforme in $\{1,\dots,i\}$; se $j\le k$ l'elemento $i$ sostituisce quello nello slot $j$, altrimenti viene scartato.

Teorema. Dopo aver letto $i\ge k$ elementi, ognuno di essi è nel campione con probabilità esattamente $k/i$ (e ogni sottoinsieme di $k$ elementi ha probabilità $1/\binom ik$).

Dimostrazione (per induzione su $i$). Per $i=k$ la probabilità è $1=k/k$. Passo $i\to i+1$: l'elemento nuovo entra con probabilità $P(j\le k)=\frac k{i+1}$. Un elemento vecchio $x$ è nel campione con probabilità $\frac ki$ (ipotesi induttiva); dato che lo è, occupa uno slot e viene espulso solo se $j$ coincide con quello slot, cioè con probabilità $\frac1{i+1}$; quindi sopravvive con probabilità $\frac i{i+1}$ e la probabilità totale è $\frac ki\cdot\frac i{i+1}=\frac k{i+1}$. $\square$ Per $k=1$: l'elemento $j$ sostituisce con probabilità $\frac1j$ e poi non viene più sostituito con probabilità $\prod_{t=j+1}^N(1-\frac1t)=\frac jN$, in totale $\frac1j\cdot\frac jN=\frac1N$.

Esempio ($N=5$, $k=2$, elenco esatto). Si enumerano tutte le scelte casuali del programma (con probabilità razionali): ognuno dei $5$ elementi è nel campione con probabilità $\frac25$ e ognuno dei $\binom52=10$ sottoinsiemi ha probabilità $\frac1{10}$. L'algoritmo "ogni nuovo elemento entra con probabilità $\frac12$ in uno slot a caso" non è uniforme: con $N=5$, $k=2$ dà $\frac{27}{64}$ ai primi due elementi e $\frac12$ all'ultimo, invece di $\frac25$ a tutti. Per $k=1$ e sostituzione con probabilità $\frac12$ i cinque elementi hanno probabilità $\frac1{16},\frac1{16},\frac18,\frac14,\frac12$.

Contatore di Morris

Problema. Contare fino a $n$ con $O(\log\log n)$ bit invece di $\log_2n$, accettando un'approssimazione.

Algoritmo (contatore di Morris). Si memorizza solo un piccolo intero $X$, inizialmente $0$. A ogni evento: con probabilità $2^{-X}$ si pone $X\leftarrow X+1$. Stima: $\tilde n=2^X-1$.

Idea: $X$ approssima $\log_2n$, e per incrementarlo da $X$ a $X+1$ servono in media $2^X$ eventi, ciascuno dei quali pesa quanto la stima.

Teorema. Dopo $n$ eventi, $E[2^{X_n}-1]=n$ (stima non distorta) e $\text{Var}(2^{X_n}-1)=\frac{n(n-1)}2$.

Dimostrazione. Sia $Z_n=2^{X_n}$. Dato $X_n=x$: $Z_{n+1}=2^{x+1}$ con probabilità $2^{-x}$ e $2^x$ altrimenti. Quindi $$E[Z_{n+1}\mid X_n=x]=2^{x+1}2^{-x}+2^x(1-2^{-x})=2^x+1,\qquad E[Z_{n+1}^2\mid X_n=x]=4^{x+1}2^{-x}+4^x(1-2^{-x})=4^x+3\cdot2^x.$$ Prendendo i valori attesi: $E[Z_{n+1}]=E[Z_n]+1$ e $E[Z_{n+1}^2]=E[Z_n^2]+3E[Z_n]$. Con $Z_0=1$ la prima dà $E[Z_n]=n+1$. Sostituendo nella seconda, $E[Z_n^2]=1+3\sum_{i=1}^ni=1+\frac{3n(n+1)}2$. Infine $\text{Var}(Z_n)=1+\frac{3n(n+1)}2-(n+1)^2=\frac{n^2-n}2$. $\square$

Esempio ($n=3$). $X$ vale $1,2,3$ con probabilità $\frac14,\frac58,\frac18$ (da $X=1$: sale con $\frac12$; da $X=2$: sale con $\frac14$). Le stime $2^X-1=1,3,7$ danno $E=\frac14+\frac{15}8+\frac78=3$, come deve. L'enumerazione esatta con frazioni conferma $E[Z_n]=n+1$ e $\text{Var}=\frac{n(n-1)}2$ per $n=1,2,3,4,10,30$.

Precisione. La deviazione standard $\approx n/\sqrt2$ è grande quanto la stima: un solo contatore è troppo grossolano. Rimedi:

• Media di $s$ contatori indipendenti: varianza $\frac{n(n-1)}{2s}$; per Chebyshev $P(|\tilde n-n|\ge\varepsilon n)\le\frac{1}{2s\varepsilon^2}$, e basta $s\ge\frac1{2\varepsilon^2\delta}$ per probabilità di errore $\delta$ (per $\varepsilon=\delta=0{,}1$: $s=500$).
• Base $1+a$: si incrementa con probabilità $(1+a)^{-X}$ e si stima $\frac{(1+a)^X-1}a$; la stima resta non distorta e $\text{Var}=\frac{a\,n(n-1)}2$ (verificato esattamente per $a=\frac12$ e $a=\frac1{10}$, $n=20$: $95$ e $19$). Con $a=2\varepsilon^2\delta$ un solo contatore dà $P(\text{errore}\ge\varepsilon n)\le\delta$ per Chebyshev. Con $\varepsilon=\delta=0{,}1$ ($a=0{,}002$) e $n=10^6$, $X\approx3800$: $12$ bit invece dei $20$ del conteggio esatto.

Simulazione ($n=1000$, $a=1$). Su $20\,000$ prove la media delle stime è $998{,}2$ e la deviazione standard $706{,}7$ (teoria $\sqrt{1000\cdot999/2}=706{,}8$); $X$ vale in media $9{,}69$ ($\log_2 1001=9{,}97$) e al massimo $13$. La media di $500$ contatori ha deviazione standard $34{,}1$ (teoria $31{,}6$), e in $300$ prove nessuna volta l'errore supera il $10\%$ di $n$.

Elementi distinti

Problema. Quanti valori diversi compaiono nel flusso? Con memoria $\Omega(d)$ è banale (insieme dei visti); si vuole una stima con memoria piccola.

Idea comune. Una funzione hash casuale manda ogni valore distinto in un punto uniforme di $(0,1)$, e i duplicati nello stesso punto: la statistica che si guarda dipende solo dall'insieme dei distinti.

Flajolet-Martin (idea)

Si mandano gli elementi in stringhe di bit e si indica con $\rho(x)$ il numero di zeri finali dell'hash: $P(\rho\ge j)=2^{-j}$. Con $d$ valori distinti, circa $d2^{-j}$ hanno $\rho\ge j$: l'indice massimo $R=\max\rho$ è vicino a $\log_2d$ e $2^R$ è una stima grossolana di $d$ (dopo una costante di correzione $\approx0{,}77$). Una sola copia è però molto variabile. In $60$ prove con $d=1000$ ($\log_2d=9{,}97$) il valore $R$ va da $8$ a $18$ (mediano $10$): per ottenere una buona stima si usano molte copie e la mediana.

Il $k$-esimo minimo

Algoritmo ($k$-esimo minimo, bottom-$k$). Con $h$ che manda i valori in $(0,1)$ si tengono i $k$ valori hash distinti più piccoli visti. Se $v_k$ è il $k$-esimo più piccolo la stima è $$\hat d=\frac{k-1}{v_k}.$$ Se i distinti sono meno di $k$ si conta direttamente.

Perché funziona. Con $d$ punti uniformi indipendenti, il minimo vale in media $\frac1{d+1}$: $P(\min>t)=(1-t)^d$, quindi $E[\min]=\int_0^1(1-t)^d\,dt=\frac1{d+1}$ (per $d=50$ la simulazione dà $0{,}01959$ contro $\frac1{51}=0{,}01961$). Più punti ci sono, più i minimi sono piccoli: $d\approx1/v_1$. Con $k=1$ la stima $1/v_1$ è però inutilizzabile: $E[1/V_1]=d\int_0^1 v^{-1}(1-v)^{d-1}dv$ diverge vicino a $v=0$, quindi media e varianza sono infinite. Per $d=50$, da $P(V_1>t)=(1-t)^{50}$, il $5\%$ e il $95\%$ delle stime valgono $1/(1-0{,}05^{1/50})=17{,}2$ e $1/(1-0{,}95^{1/50})=975$.

Teorema. Per $k\ge3$, $E\bigl[\frac{k-1}{V_k}\bigr]=d$ (stima non distorta) e lo scarto relativo è circa $\frac1{\sqrt{k-2}}$.

Dimostrazione della media. La densità del $k$-esimo più piccolo di $d$ punti uniformi è $\frac{d!}{(k-1)!(d-k)!}v^{k-1}(1-v)^{d-k}$ (si sceglie quali $k-1$ punti stanno sotto, quale è $V_k$, e gli altri $d-k$ stanno sopra). Quindi $$E\Bigl[\frac1{V_k}\Bigr]=\frac{d!}{(k-1)!(d-k)!}\int_0^1v^{k-2}(1-v)^{d-k}dv=\frac{d!}{(k-1)!(d-k)!}\cdot\frac{(k-2)!(d-k)!}{(d-1)!}=\frac d{k-1},$$ dove l'integrale è la funzione Beta. $\square$ La varianza si calcola allo stesso modo e dà $\text{sd}(\hat d)\approx d/\sqrt{k-2}$. Per Chebyshev, con $k=\lceil\frac{2}{\varepsilon^2}\rceil+2$ la stima è entro $(1\pm\varepsilon)d$ con probabilità $\ge\frac12$, amplificabile con la mediana di copie indipendenti. La memoria è $O(k)$ valori hash, indipendente da $d$ e da $N$.

Simulazione ($d=10^5$, flusso con ripetizioni). Con $k=100$ su $200$ prove: media $99\,852$, scarto relativo $0{,}104$ (teoria $0{,}101$), $95\%$ delle stime entro il $20\%$. Con $k=400$: media $100\,107$, scarto $0{,}051$ (teoria $0{,}050$), $96\%$ entro il $10\%$.

Count-Min sketch

Problema. Stimare la frequenza $f(x)$ di ogni elemento $x$ (numero di arrivi) con memoria piccola, per esempio per trovare gli elementi più frequenti.

Definizione (Count-Min sketch). Una matrice $C$ di $d$ righe e $w$ colonne di contatori, inizialmente $0$, e $d$ funzioni hash indipendenti $h_1,\dots,h_d:U\to[w]$ da una famiglia universale. Arrivo di $x$: $C[j][h_j(x)]\mathrel{+}=1$ per ogni riga $j$. Stima: $\hat f(x)=\min_{j}C[j][h_j(x)]$.

Teorema (garanzia $(\varepsilon,\delta)$). Sia $F=\sum_yf(y)=N$ il numero totale di arrivi. Con $w=\lceil e/\varepsilon\rceil$ e $d=\lceil\ln(1/\delta)\rceil$: $$f(x)\le\hat f(x)\quad\text{sempre},\qquad P\bigl(\hat f(x)\le f(x)+\varepsilon N\bigr)\ge1-\delta.$$

Dimostrazione. Nella riga $j$ la cella di $x$ contiene $f(x)$ più l'errore $Y_j=\sum_{y\ne x}f(y)\,\mathbf 1[h_j(y)=h_j(x)]\ge0$, quindi $\hat f(x)=f(x)+\min_jY_j\ge f(x)$. Per l'universalità $P(h_j(y)=h_j(x))\le\frac1w$ e, per la linearità, $E[Y_j]\le\frac Nw\le\frac{\varepsilon N}e$. Markov (in modo valido perché $Y_j\ge0$): $P(Y_j\ge\varepsilon N)\le\frac1e$. Le righe sono indipendenti, e $\hat f(x)>f(x)+\varepsilon N$ richiede che tutte le righe superino $\varepsilon N$: probabilità $\le e^{-d}\le\delta$. $\square$

Esempio piccolo ($p=17$, $w=5$, $d=3$). Funzioni $h_j(x)=((a_jx+b_j)\bmod17)\bmod5$ con $(a_j,b_j)=(3,4),(5,1),(7,9)$; flusso di $N=20$ arrivi con frequenze $f(3)=7$, $f(8)=5$, $f(12)=3$, $f(15)=2$, $f(9)=2$, $f(4)=1$. Posizioni (riga $0$/$1$/$2$): $3\mapsto3/1/3$, $8\mapsto1/2/4$, $12\mapsto1/0/3$, $15\mapsto0/3/2$, $9\mapsto4/2/4$, $4\mapsto1/4/3$. Contatori finali:

La colonna $1$ della riga $0$ vale $5+3+1=9$ (collidono $8,12,4$), la colonna $2$ della riga $1$ vale $5+2=7$ ($8$ e $9$), la colonna $3$ della riga $2$ vale $7+3+1=11$. Stime: $\hat f(8)=\min(9,7,7)=7$ (vero $5$: sovrastima di $2$, perché $8$ collide con $9$ nelle righe $1$ e $2$ e con $12$ e $4$ nella riga $0$), tutte le altre sono esatte; l'elemento mai visto $6$ va nelle colonne $0,4,0$ ($3\cdot6+4=22\equiv5$, $5\cdot6+1=31\equiv14$, $7\cdot6+9=51\equiv0$ modulo $17$, poi modulo $5$), con contatori $2,1,0$, e ha stima $\min(2,1,0)=0$. Nessuna stima è sotto la frequenza vera.

Dimensionamento. Per $\varepsilon=0{,}001$ e $\delta=0{,}01$: $w=\lceil e/0{,}001\rceil=2719$ colonne e $d=\lceil\ln100\rceil=5$ righe, cioè $13\,595$ contatori indipendentemente dal numero di elementi distinti. L'errore additivo è $\le\varepsilon N$: va bene per gli elementi frequenti ($f(x)\gg\varepsilon N$), non per quelli rari.

Simulazione. Su un flusso di $N=200\,000$ arrivi su $20\,000$ valori con distribuzione di tipo Zipf, con $\varepsilon=0{,}01$, $\delta=0{,}05$ ($w=272$, $d=3$) l'errore massimo osservato su tutti i valori è $1357<\varepsilon N=2000$; con $\varepsilon=0{,}001$, $\delta=0{,}01$ ($w=2719$, $d=5$) è $44<200$.

import random, heapq, hashlib
from math import e, log, ceil

def riserva(flusso, k):
    """Campione uniforme di k elementi da un flusso di lunghezza ignota."""
    r = []
    for i, x in enumerate(flusso, 1):
        if i <= k:
            r.append(x)
        elif (j := random.randrange(1, i + 1)) <= k:   # entra con probabilità k/i
            r[j - 1] = x                               # al posto di uno slot a caso
    return r

def morris(n):
    X = 0
    for _ in range(n):
        if random.random() < 2.0 ** -X:                # incrementa con probabilità 2^-X
            X += 1
    return 2 ** X - 1                                  # stima non distorta di n

def h01(x):                                            # hash in (0,1)
    d = hashlib.blake2b(str(x).encode(), digest_size=8).digest()
    return (int.from_bytes(d, "big") + 0.5) / 2 ** 64

def distinti(flusso, k):
    """k-esimo minimo: stima (k-1)/v_k del numero di elementi distinti."""
    heap, visti = [], set()                            # max-heap (valori negati) dei k minimi
    for x in flusso:
        v = h01(x)
        if v in visti:
            continue
        if len(heap) < k:
            heapq.heappush(heap, -v); visti.add(v)
        elif v < -heap[0]:
            visti.discard(-heapq.heappushpop(heap, -v)); visti.add(v)
    return len(heap) if len(heap) < k else (k - 1) / -heap[0]

class CountMin:
    P = (1 << 61) - 1
    def __init__(self, eps, delta):
        self.w, self.d = ceil(e / eps), ceil(log(1 / delta))
        self.ab = [(random.randrange(1, self.P), random.randrange(self.P)) for _ in range(self.d)]
        self.C = [[0] * self.w for _ in range(self.d)]
    def _col(self, j, x):
        a, b = self.ab[j]
        return ((a * x + b) % self.P) % self.w
    def aggiungi(self, x, c=1):
        for j in range(self.d):
            self.C[j][self._col(j, x)] += c
    def stima(self, x):                                # mai sotto la frequenza vera
        return min(self.C[j][self._col(j, x)] for j in range(self.d))

Con $200\,000$ arrivi su $50\,000$ valori casuali ($49\,080$ distinti) il programma completo stima $48\,401$ elementi distinti con $k=400$, e un valore comparso $4$ volte ha stima $73$ in un Count-Min con $\varepsilon=0{,}001$ (errore $69<\varepsilon N=200$: per gli elementi rari la garanzia additiva è debole).

Domande d'esame

1. Dimostra che il reservoir sampling con $k=1$ restituisce un elemento uniforme. L'elemento $j$ entra con probabilità $\frac1j$ e sopravvive ai passi $j+1,\dots,N$ con probabilità $\prod_{t=j+1}^N\frac{t-1}t=\frac jN$; prodotto $\frac1N$.
2. Dimostra che il contatore di Morris è non distorto. $E[2^{X_{n+1}}\mid X_n]=2^{X_n}+1$, quindi $E[2^{X_n}]=1+n$ e $E[2^{X_n}-1]=n$.
3. Con quanti contatori di Morris si ottiene errore relativo $\le10\%$ con probabilità $\ge90\%$? Chebyshev: $\frac{n(n-1)/(2s)}{(0{,}1n)^2}<\frac1{0{,}02s}\le0{,}1$, quindi $s\ge500$.
4. Perché nello sketch Count-Min si prende il minimo e non la media delle righe? Ogni cella contiene la frequenza vera più un errore $\ge0$ da collisioni: il minimo è la cella meno contaminata e non è mai sotto $f(x)$; la media includerebbe anche le righe peggiori.