Revisiting Bayesian Variable Selection via Optimization

Resumen Técnico: Revisión de la Selección de Variables Bayesiana mediante Optimización

1. Planteamiento del Problema

La selección de variables en regresión lineal es un tema central en la investigación estadística, especialmente cuando el número de predictores ($p$) es comparable o superior al tamaño de la muestra ($n$). Los métodos bayesianos tradicionales, que utilizan priores de contracción (como spike-and-slab o priores globales-locales como el horseshoe), han tenido gran éxito mediante el uso de Muestreo de Gibbs colapsado para integrar coeficientes y varianza del ruido.

Sin embargo, persiste una enigmas fundamental: la verosimilitud marginal (obtenida integrando los coeficientes de regresión $\theta$ y la varianza del ruido $\sigma^2$) no es log-cóncava. Esto implica que:

• No existe garantía teórica de que los algoritmos de optimización encuentren el óptimo global.
• La búsqueda del modo de la distribución posterior (estimación de máxima verosimilitud marginal o tipo-II) podría quedar atrapada en óptimos locales.
• Existe incertidumbre sobre la fiabilidad de los métodos existentes para encontrar el mejor modelo de selección de variables.

El objetivo del artículo es abordar este problema desde la perspectiva de la optimización no convexa, demostrando que, a pesar de la falta de log-concavidad, el problema es mucho más benigno de lo que parece y puede resolverse globalmente y eficientemente.

2. Metodología

El autor propone un enfoque basado en la optimización de la verosimilitud marginal (enfoque de máxima verosimilitud tipo-II) utilizando una estructura matemática específica.

• Formulación del Problema:
  Se define una función de pérdida como el logaritmo negativo de la verosimilitud marginal $L(y|D)$, donde $D = \text{diag}(d_1, \dots, d_p)$ son los parámetros de precisión que controlan la contracción de cada coeficiente $\theta_j$.
  $$\text{Loss}(D) = -\log L(y|D)$$
  El problema se formula como la minimización de esta función sobre un conjunto convexo y compacto de $D$.

• Estructura Diferencia de Funciones Convexas (DC):
  La contribución teórica central es demostrar que la función de pérdida puede reescribirse como la diferencia de dos funciones convexas ($g(D) - h(D)$):

  • $g(D)$: Una función convexa que incluye términos logarítmicos de los parámetros de precisión y restricciones de dominio.
  • $h(D)$: Una función convexa derivada del término de ajuste de datos y la varianza.

• Algoritmo Propuesto (Algoritmo DC):
  Se emplea un algoritmo iterativo simple basado en la aproximación de primer orden de la parte cóncava ($h$):

  1. Inicializar $D^{(0)}$.
  2. En cada iteración $t$, aproximar $h(D)$ linealmente alrededor de $D^{(t)}$.
  3. Resolver el subproblema convexo resultante para obtener $D^{(t+1)}$.
  4. Para problemas de alta dimensión, se utiliza la fórmula de Woodbury para reducir la complejidad computacional de la inversión de matrices a $O(\min(n, p)^3)$.
  5. Se incluye una estrategia de exploración de fronteras para acelerar la convergencia hacia los límites del dominio (donde ocurre la selección de variables, es decir, $d_j \to \infty$).

• Extensiones:
  El marco se extiende fácilmente a esparsidad estructurada (grupos, superposiciones, ordenamiento radial) imponiendo restricciones adicionales en el conjunto convexo $D$, resolviéndose mediante proyecciones alternas o regresión isotónica.

3. Contribuciones Clave

1. Garantía de Optimalidad Global: Bajo condiciones suaves (conjunto compacto y cerrado), el algoritmo DC converge al óptimo global a una tasa lineal. Esto resuelve la duda sobre la fiabilidad de encontrar el modo de la verosimilitud marginal en selección de variables bayesiana.
2. Análisis de Convergencia: Se demuestra teóricamente que los puntos estacionarios no óptimos (que podrían ser trampas para otros algoritmos) son evitables o están en la frontera, y el algoritmo los supera.
3. Eficiencia y Simplicidad: El algoritmo es:
   • Libre de ajuste (tuning-free): No requiere selección de hiperparámetros complejos.
   • Rápido: Convierte un problema no convexo en una serie de subproblemas convexos fáciles de resolver.
   • Escalable: Funciona bien en dimensiones altas ($p \gg n$).
4. Alternativa al MCMC: Se propone como una alternativa eficiente o como un "punto de inicio" (warm-start) para algoritmos MCMC tradicionales, reduciendo significativamente el tiempo de burn-in.

4. Resultados

El artículo valida el método mediante estudios numéricos y aplicaciones reales:

• Simulaciones Sintéticas:

  • Se comparó el algoritmo DC contra el algoritmo EM para Determinación de Relevancia Automática (ARD) y el Descenso de Gradiente Proyectado (PGD).
  • Precisión: DC logró una selección de variables perfecta (tasa de verdaderos positivos = 1, falsos positivos = 0) en todos los escenarios, confirmando la garantía de optimalidad global.
  • Velocidad: DC requirió el menor número de iteraciones (aprox. 80-300) y fue el más rápido en tiempo de ejecución (menos de 2.3 segundos por réplica en CPU estándar), superando a ARD y PGD.

• Datos de Alto Dimensionalidad (Dorothea):

  • Aplicado a un conjunto de datos de descubrimiento de fármacos con $n=1,150$ y $p=91,598$.
  • El algoritmo convergió rápidamente, eliminando el 99.97% de las características (empujando sus $d_j$ a la frontera superior), seleccionando un número de características comparable a Lasso y Group Lasso, pero incorporando la incertidumbre de la varianza del ruido.

• Aplicación en Datos Espaciales (Terremotos Ridgecrest):

  • Se utilizó para cuantificar el riesgo de réplicas tras el terremoto de 2019 en California.
  • Se impuso una contracción monótona radial (mayor contracción a mayor distancia de la falla).
  • El algoritmo convergió en solo 35 iteraciones (0.20 segundos), identificando patrones espaciales coherentes con la física del fenómeno (actividad concentrada en la falla, disminuyendo con la distancia).

5. Significado e Impacto

Este trabajo es significativo porque:

• Desmitifica la optimización en Bayesiana: Demuestra que la falta de log-concavidad en la verosimilitud marginal no es un obstáculo insuperable para encontrar soluciones globales, cambiando la percepción sobre la dificultad de la selección de variables bayesiana.
• Puente entre Optimización y Estadística Bayesiana: Ofrece un método determinista, rápido y teóricamente garantizado que complementa o reemplaza a los métodos estocásticos (MCMC) en muchos contextos, especialmente en alta dimensión.
• Versatilidad: La capacidad de incorporar fácilmente restricciones de esparsidad estructurada lo hace aplicable a problemas complejos de aprendizaje automático y análisis espacial.
• Reproducibilidad: El código fuente es público, facilitando la adopción por parte de la comunidad estadística y de ciencia de datos.

En conclusión, el artículo presenta un avance teórico y práctico que transforma un problema de optimización no convexa aparentemente difícil en un procedimiento eficiente y robusto para la selección de variables bayesiana.