Revisiting Bayesian Variable Selection via Optimization

Probleemstelling

Variabele selectie in lineaire regressie is een kernonderwerp in de statistiek, vooral wanneer het aantal potentiële voorspellers ($p$) vergelijkbaar is met of groter is dan het steekproefgrootte ($n$). Bayesiaanse methoden zijn populair omdat ze onzekerheid in zowel regressiecoëfficiënten als ruisvariatie kunnen modelleren, vaak via schuine priors (zoals spike-and-slab of global-local shrinkage priors) en efficiënte MCMC-algoritmen (zoals collapsed Gibbs sampling).

Echter, er blijft een fundamenteel probleem bestaan: de marginale likelihood (verkregen door integratie van de regressiecoëfficiënten en de ruisvariatie) is niet log-concaaf. Dit betekent dat er geen theoretische garantie is dat optimalisatie-algoritmen het globale optimum vinden; ze kunnen vastlopen in lokale optima. Dit creëert twijfel over de betrouwbaarheid van schattingen zoals de Type-II maximum likelihood (empirische Bayes) of de modus van de posterior voor de schuine parameters.

Methodologie

De auteur benadert dit probleem vanuit het perspectief van niet-convexe optimalisatie. De kern van de methode is als volgt:

1. Formulering als Verliesfunctie:
   De negatieve log-marginale likelihood wordt behandeld als een verliesfunctie voor de latente precisieparameters $D = \text{diag}(d_1, \dots, d_p)$.
   $$\mathcal{L}(D) = -\frac{1}{2}\log|D| + \frac{1}{2}\log|X^\top X + D| + (a_0 + \frac{n}{2}) \log \left( b_0 + \frac{1}{2} \text{RSS}(D) \right)$$
   Waarbij $\text{RSS}(D)$ de residu som van kwadraten is na integratie van $\theta$ en $\sigma^2$.

2. Difference-of-Convex (DC) Structuur:
   De verliesfunctie wordt herschreven als het verschil van twee convexe functies: $f(D) = g(D) - h(D)$.

   • $g(D)$ is convex (bevat $-\log|D|$).
   • $h(D)$ is convex (bevat termen die voortvloeien uit de inverse van de covariance matrix en de log-term).
     Deze structuur maakt het mogelijk om het DC-algoritme (ook wel DCA genoemd) toe te passen.

3. Het DC-Algoritme:
   Het algoritme is een iteratief procedure:

   • Start met een initiële $D^{(0)}$.
   • Benader $h(D)$ lineair rond het huidige punt $D^{(t)}$ (via een eerste-orde Taylor-ontwikkeling).
   • Los het onderliggende convexe optimalisatieprobleem op om $D^{(t+1)}$ te vinden.
   • Herhaal tot convergentie.
     De update is eenvoudig en vereist geen tuning van hyperparameters (zoals leerstappen). Voor box-beperkingen ($a \le d_j \le b$) heeft de update een gesloten vorm.

4. Uitbreiding tot Gestructureerde Sparsiteit:
   De methode is flexibel genoeg om extra convexe beperkingen op te leggen aan $D$ om gestructureerde sparsiteit te induceren (bijv. groepen, overlappende groepen, of radiale monotoniciteit). Dit wordt gedaan via projectie-algoritmen (zoals isotone regressie) binnen elke iteratie.

Belangrijkste Bijdragen

• Theoretische Garantie: De auteur bewijst dat onder milde voorwaarden (compacte verzameling, lokale kwadratische groei) het DC-algoritme convergeert naar het globale optimum met een lineaire snelheid. Dit weerlegt de angst dat het gebrek aan log-concaafheid noodzakelijk leidt tot lokale optima; het probleem is "minder kwaadaardig" dan gedacht.
• Efficiëntie en Eenvoud: Het algoritme is eenvoudig te implementeren, vereist geen tuning en is schaalbaar. Het kan dienen als een efficiënt alternatief of als een "warm-start" voor traditionele MCMC-methoden.
• Uitbreiding naar Posterior Modus: De theorie geldt niet alleen voor Type-II maximum likelihood, maar ook voor het vinden van de modus van de posterior onder bepaalde priors.
• Open Source: De code is publiek beschikbaar, wat reproduceerbaarheid bevordert.

Resultaten

De methode werd getest op synthetische data, een hoogdimensionale dataset (Dorothea) en een ruimtelijke toepassing.

1. Synthetische Data:

   • Vergelijking met EM-algoritmen (ARD) en Projected Gradient Descent (PGD).
   • Het DC-algoritme bereikte in alle scenario's (tot $p=5000$) een perfecte variabele selectie (TPR=1, FPR=0), wat de theoretische garantie van globale optimaliteit bevestigt.
   • Het vereiste het minste aantal iteraties en was het snelst in rekentijd (onder 2,3 seconden per replicatie op een Macbook Pro).

2. Dorothea Dataset (Hoogdimensionaal):

   • Toepassing op $n=1.150$ en $p=91.598$ (drug discovery).
   • Na 5.000 iteraties (op GPU) werden 99,97% van de features naar de bovenste grens geduwd (effectief geselecteerd als ruis).
   • De methode selecteerde een vergelijkbaar aantal relevante features als Lasso en Group Lasso, maar doet dit met een volledig Bayesiaanse behandeling van onzekerheid.

3. Ruimtelijke Toepassing (Aardschokken):

   • Modellering van naschokken na de Ridgecrest-aardschokken (2019).
   • Toepassing van radiale monotone schuining: parameters verder van de breuklijn worden sterker naar nul geschaald.
   • Het algoritme convergeerde in slechts 35 iteraties (0,20 seconden) en produceerde een gladde, fysiek interpretabele schatting van de activiteit die overeenkomt met de verwachte afname van activiteit met afstand tot de breuk.

Betekenis en Conclusie

Dit artikel biedt een doorbraak in het begrip van Bayesiaanse variabele selectie. Het toont aan dat het optimalisatieprobleem achter de marginale likelihood, ondanks de niet-convexe aard, goed hanteerbaar is via DC-optimalisatie.

De belangrijkste implicaties zijn:

• Betrouwbaarheid: Het globale optimum kan betrouwbaar worden gevonden zonder afhankelijk te zijn van MCMC-kettingen die mogelijk niet convergeren of lange brand-in periodes vereisen.
• Snelheid: Het algoritme is extreem snel en schaalbaar naar zeer hoge dimensies ($p \gg n$).
• Flexibiliteit: Het biedt een natuurlijke weg om complexe structurele beperkingen (zoals ruimtelijke of groepsgewijze afhankelijkheid) toe te passen binnen een Bayesiaans raamwerk.

De methode positioneert zich als een krachtig alternatief voor MCMC, vooral geschikt als standalone oplossing voor modus-vindingsproblemen of als snelle initialisatie voor verdere Bayesiaanse inferentie.