Mínimos Cuadrados No Negativos

Imagine que uno tiene una matriz de datos \(X \in \mathbb{R}^{n \times p}\) que consiste en \(n\) observaciones, cada una con \(p\) características, así como un vector de respuesta \(y \in \mathbb{R}^n\). Queremos construir un modelo para \(y\) usando las columnas de características en \(X\). En los mínimos cuadrados ordinarios (MCO), uno busca un vector de coeficientes \(\hat{\beta} \in \mathbb{R}^p\) tal que \[\begin{aligned} \hat{\beta} = \underset{\beta \in \mathbb{R}^p}{\text{argmin}} \quad \|y - X\beta \|_2^2. \end{aligned}\]

En los mínimos cuadrados no negativos (Non Negative Least Squares, NNLS), buscamos coeficientes \(\hat{\beta} \in \mathbb{R}^p\) vectoriales de manera que se minimice \(\quad \|y - X\beta \|_2^2\) sujeto al requisito adicional de que cada elemento de no sea negativo. Hay varias formas de realizar NNLS en R. Se presentan los dos usuales con los paquetes nnls y glmnet, que suelen tener problemas, y al final el de Kjytay (2019), quien propone un proceso sencillo y que corre sin mayores problema al implementarlo. Generemos algunos datos falsos que usaremos para el resto de la publicación:

set.seed(1)
n <- 100; p <- 10
x <- matrix(rnorm(n * p), nrow = n)
y <- x %*% matrix(rep(c(1, -1), length.out = p), ncol = 1) + rnorm(n)

Método 1: el paquete nnls

library(nnls)
mod1 <- nnls(x, y)
mod1$x 

##  [1] 0.9073423 0.0000000 1.2971069 0.0000000 0.9708051 0.0000000 1.2002310
##  [8] 0.0000000 0.3947028 0.0000000

Método 2: el paquete glmnet

La función glmnet() resuelve el problema de minimización \[\begin{aligned} \hat{\beta} = \underset{\beta \in \mathbb{R}^p}{\text{argmin}} \quad \frac{1}{2n} \|y - X\beta\|_2^2 + \lambda \left[\frac{1-\alpha}{2}\|\beta\|_2^2 + \alpha \|\beta\|_1 \right], \end{aligned}\] donde \(\alpha\) y \(\lambda\) son hiperparámetros que elige el usuario. Al fijar \(\alpha = 1\) (el valor predeterminado) y \(\lambda = 0\), glmnet() termina resolviendo el problema de los Mínimos Cuadrados Ordinarios (Ordinary Least Sqaures, OLS). Al establecer lower.limits = 0, esto obliga a que los coeficientes sean no negativos. También debemos establecer intercept = FALSE para no tener un término de intercepción extraño.

library(glmnet)

## Loading required package: Matrix

## Loaded glmnet 3.0-1

mod2 <- glmnet(x, y, lambda = 0, lower.limits = 0, intercept = FALSE)
coef(mod2)

## 11 x 1 sparse Matrix of class "dgCMatrix"
##                    s0
## (Intercept) .        
## V1          0.9073427
## V2          .        
## V3          1.2971070
## V4          .        
## V5          0.9708049
## V6          .        
## V7          1.2002310
## V8          .        
## V9          0.3947028
## V10         .

Método 3: el paquete bvls

Los NNLS son un caso especial de Mínimos Cuadrados de Variable Acotada (Bounded Variable Least Squares, BVLS), donde en lugar de tener restricciones \(\beta_j\geq0\) para cada \(j = 1, 2, \dots, p\), uno tiene restricciones \(a_j\leq\beta_j\leq b_j\) para cada \(j\). Los BVLS se implementan en la función bvls() del paquete bvls:

library(bvls)
mod3 <- bvls(x, y, bl = rep(0, p), bu = rep(Inf, p))
mod3$x

##  [1] 0.9073423 0.0000000 1.2971069 0.0000000 0.9708051 0.0000000 1.2002310
##  [8] 0.0000000 0.3947028 0.0000000

En lo anterior, bl contiene los límites inferiores para los coeficientes, mientras que bu contiene los límites superiores.

Referencias

Cosas que pensé en un punto. Tres formas de hacer mínimos cuadrados no negativos en R.