Revisiting Bayesian Variable Selection via Optimization

1. Il Problema

La selezione delle variabili nella regressione lineare è un tema centrale nella ricerca statistica, specialmente quando il numero di predittori ($p$) è comparabile o superiore alla dimensione del campione ($n$). I metodi bayesiani sono ampiamente utilizzati per gestire l'incertezza sia sui coefficienti di regressione che sulla varianza del rumore, impiegando prior di shrinkage (discrete o continue) e campionatori di Gibbs collassati.

Tuttavia, persiste un enigma fondamentale: la marginal likelihood (ottenuta integrando fuori i coefficienti $\theta$ e la varianza $\sigma^2$) non è log-concava. Di conseguenza, non vi è alcuna garanzia teorica che gli algoritmi di ottimizzazione standard (come la massimizzazione della verosimiglianza marginale o la ricerca del modo della posterior) riescano a trovare l'ottimo globale, rischiando di rimanere intrappolati in ottimi locali. Questo dubbio sulla affidabilità della ricerca del modo globale ha limitato la percezione della robustezza dei metodi di selezione bayesiana.

2. Metodologia

L'autore affronta il problema da una prospettiva di ottimizzazione non convessa, proponendo un approccio basato sulla struttura Differenza di Funzioni Convesse (DC - Difference of Convex functions).

• Formulazione del Problema:
  Considerando il modello gerarchico standard con prior normale e varianza inversa-gamma, l'autore definisce la funzione di perdita come il logaritmo negativo della marginal likelihood rispetto ai parametri di precisione $D = \text{diag}(d_1, \dots, d_p)$.
  La funzione di perdita $f(D)$ viene riscritta come:
  $$f(D) = g(D) - h(D)$$
  dove:

  • $g(D) = -\frac{1}{2}\sum \log d_j$ è una funzione convessa.
  • $h(D)$ è una funzione che include termini logaritmici e quadratici derivanti dalla likelihood, che si dimostra essere convessa (o la sua parte negativa è concava).

• L'Algoritmo DC:
  Viene proposto un algoritmo iterativo semplice e privo di parametri di tuning (tuning-free):

  1. Si parte da un punto iniziale $D^{(0)}$.
  2. All'iterazione $t$, la funzione convessa $h(D)$ viene approssimata con il suo sviluppo di Taylor del primo ordine (approssimazione affine) attorno a $D^{(t)}$.
  3. Si risolve un sottoproblema di ottimizzazione convessa per aggiornare $D$:
     $$D^{(t+1)} = \arg\min_{D \in \mathcal{D}} \left\{ g(D) - \left( h(D^{(t)}) + \langle \nabla h(D^{(t)}), D - D^{(t)} \rangle \right) \right\}$$
  4. Per vincoli di scatola (box constraints) $a \le d_j \le b$, l'aggiornamento ammette una forma chiusa che coinvolge la proiezione sul dominio e il gradiente di $h$.

• Gestione della Sparsità Strutturata:
  Il metodo è estendibile a vincoli di sparsità strutturata (gruppi, sovrapposizioni, monotonia radiale) utilizzando algoritmi di proiezione alternata o regressione isotona generalizzata (es. algoritmo PAVA).

3. Contributi Chiave e Risultati Teorici

• Convergenza Globale:
  Il contributo teorico principale è la dimostrazione che, sotto condizioni di compattezza e convessità del dominio, l'algoritmo DC converge all'ottimo globale a un tasso lineare.

  • Viene stabilito che i punti stazionari non ottimali (che potrebbero intrappolare algoritmi meno robusti) sono esclusi o si trovano sul bordo del dominio, dove il gradiente non è nullo.
  • Viene introdotto un teorema che garantisce la convergenza globale assumendo condizioni di crescita quadratica locale e bound sull'errore residuo (condizioni A1, A2, A3 nel paper).

• Superiorità rispetto agli Ottimi Locali:
  A differenza della massimizzazione coordinata (coordinate-wise maximization) che non garantisce l'ottimalità globale in problemi non convessi, l'algoritmo DC proposto trova l'ottimo globale indipendentemente dal valore iniziale, rendendo il problema della ricerca del modo ("mode finding") molto più benigno di quanto suggerito dalla mancanza di log-concavità.

• Estensione alla Posterior:
  I risultati si estendono alla massimizzazione della posterior marginale (Maximum Marginal Posterior) quando le prior aggiuntive su $D$ possono essere espresse come differenza di funzioni convesse.

4. Risultati Sperimentali

• Simulazioni:

  • Confronto con l'algoritmo EM per la Determinazione Automatica della Rilevanza (ARD) e la Discesa del Gradiente Proiettata (PGD).
  • L'algoritmo DC ha mostrato la convergenza più rapida in termini di iterazioni e tempo di calcolo (sotto 2.3 secondi per replica su CPU standard), con varianza quasi nulla tra le repliche.
  • In tutti gli scenari (da $p=500$ a $p=5000$), DC ha raggiunto una selezione perfetta delle variabili (TPR=1, FPR=0), confermando la garanzia teorica di ottimalità globale.

• Dati Reali (Dorothea):

  • Applicazione su un dataset di scoperta di farmaci con $n=1,150$ e $p=91,598$.
  • L'algoritmo ha identificato 26 caratteristiche rilevanti, un numero comparabile a Lasso e Group Lasso, ma incorporando l'incertezza su $\sigma^2$ e $\theta$.
  • La funzione di perdita è diminuita monotonamente e rapidamente.

• Applicazione Spaziale (Terremoti di Ridgecrest):

  • Modellazione del rischio di scosse di assestamento post-terremoto (2019) con vincoli di monotonia radiale (maggiore shrinkage all'aumentare della distanza dalla faglia).
  • L'algoritmo ha converguto in soli 35 iterazioni (0.20 secondi), producendo una mappa di attività sismica che rispetta le aspettative fisiche di decadimento con la distanza.

5. Significato e Implicazioni

1. Alternativa Efficiente all'MCMC:
   L'algoritmo offre un'alternativa rapida e deterministica ai metodi MCMC (come il Gibbs sampling collassato) per la selezione delle variabili. Può essere utilizzato come soluzione finale o come "warm-start" eccellente per ridurre il periodo di burn-in negli algoritmi MCMC, specialmente in problemi ad alta dimensionalità.

2. Risoluzione di un Enigma Teorico:
   Il lavoro chiarisce che, nonostante la non-concavità della marginal likelihood, la struttura specifica del problema permette una convergenza globale affidabile, rassicurando sulla solidità teorica dei metodi di selezione bayesiana basati sulla massimizzazione.

3. Flessibilità e Scalabilità:
   La metodologia è priva di parametri di tuning, facile da implementare e scalabile a problemi con $p \gg n$ (grazie all'uso della formula di Woodbury e accelerazioni GPU). Inoltre, la capacità di gestire vincoli di sparsità strutturata la rende adatta a problemi complessi come quelli spaziali o di dati con struttura a gruppi.

In sintesi, il paper trasforma un problema di ottimizzazione non convessa apparentemente intrattabile in un processo efficiente e garantito, fornendo uno strumento pratico e teorico solido per la selezione bayesiana delle variabili.