El lazo se ajusta a través de LARS (un proceso iterativo, que comienza en alguna estimación inicial $ \beta ^0$ ). Por defecto $ \beta ^0=0_p$ pero puedes cambiar esto en la mayoría de las implementaciones (y reemplazarlo por el óptimo $ \beta ^*_{old}$ ya lo has hecho). El más cercano $ \beta ^*_{old}$ es $ \beta_ {new}^*$ cuanto menor sea el número de iteración de LARS que tendrás que pisar para llegar a $ \beta_ {new}^*$ .
EDITAR:
Debido a los comentarios de user2763361 Añado más detalles a mi respuesta original.
De los comentarios que aparecen a continuación deduzco que el usuario2763361 sugiere complementar mi respuesta original para convertirla en una que pueda ser utilizada directamente (fuera de los estantes) y que al mismo tiempo sea muy eficiente.
Para hacer la primera parte, ilustraré la solución que propongo paso a paso en un ejemplo de juguete. Para satisfacer la segunda parte, lo haré usando un solucionador de puntos interiores reciente y de alta calidad. Esto se debe a que es más fácil obtener una implementación de alto rendimiento de la solución que propongo utilizando una biblioteca que puede resolver el problema del lazo por el enfoque del punto interior en lugar de tratar de hackear el algoritmo LARS o simplex para iniciar la optimización desde un punto de partida no estándar (aunque ese segundo lugar también es posible).
Nótese que a veces se afirma (en libros antiguos) que el enfoque de punto interior para resolver programas lineales es más lento que el enfoque simplex y eso puede haber sido cierto hace mucho tiempo, pero generalmente no es cierto hoy en día y ciertamente no es cierto para problemas de gran escala (esta es la razón por la que la mayoría de las bibliotecas profesionales como cplex utilizan el algoritmo de punto interior) y la pregunta es, al menos implícitamente, sobre los problemas a gran escala. Nótese también que el solucionador de puntos interiores que utilizo maneja completamente las matrices dispersas, por lo que no creo que haya una gran diferencia de rendimiento con LARS (una motivación original para utilizar LARS fue que muchos solucionadores de LP populares en ese momento no manejaban bien las matrices dispersas y estos son un rasgo característico del problema de LASSO).
Una (muy) buena implementación de código abierto del algoritmo de punto interior es ipopt en el COIN-OR biblioteca. Otra razón por la que usaré ipopt es que tiene una interfaz R, ipoptr . Encontrará una guía de instalación más exhaustiva aquí abajo doy los comandos estándar para instalarlo en ubuntu .
en el bash hacer:
sudo apt-get install gcc g++ gfortran subversion patch wget
svn co https://projects.coin-or.org/svn/Ipopt/stable/3.11 CoinIpopt
cd ~/CoinIpopt
./configure
make
make install
Entonces, como raíz, en R hacer (asumo que svn ha copiado el archivo de subversión en ~/ como lo hace por defecto):
install.packages("~/CoinIpopt/Ipopt/contrib/RInterface",repos=NULL,type="source")
A partir de aquí, estoy dando un pequeño ejemplo (mayormente del ejemplo del juguete dado por Jelmer Ypma como parte de su R envoltura a ipopt ):
library('ipoptr')
# Experiment parameters.
lambda <- 1 # Level of L1 regularization.
n <- 100 # Number of training examples.
e <- 1 # Std. dev. in noise of outputs.
beta <- c( 0, 0, 2, -4, 0, 0, -1, 3 ) # "True" regression coefficients.
# Set the random number generator seed.
ranseed <- 7
set.seed( ranseed )
# CREATE DATA SET.
# Generate the input vectors from the standard normal, and generate the
# responses from the regression with some additional noise. The variable
# "beta" is the set of true regression coefficients.
m <- length(beta) # Number of features.
A <- matrix( rnorm(n*m), nrow=n, ncol=m ) # The n x m matrix of examples.
noise <- rnorm(n, sd=e) # Noise in outputs.
y <- A %*% beta + noise # The outputs.
# DEFINE LASSO FUNCTIONS
# m, lambda, y, A are all defined in the ipoptr_environment
eval_f <- function(x) {
# separate x in two parts
w <- x[ 1:m ] # parameters
u <- x[ (m+1):(2*m) ]
return( sum( (y - A %*% w)^2 )/2 + lambda*sum(u) )
}
# ------------------------------------------------------------------
eval_grad_f <- function(x) {
w <- x[ 1:m ]
return( c( -t(A) %*% (y - A %*% w),
rep(lambda,m) ) )
}
# ------------------------------------------------------------------
eval_g <- function(x) {
# separate x in two parts
w <- x[ 1:m ] # parameters
u <- x[ (m+1):(2*m) ]
return( c( w + u, u - w ) )
}
eval_jac_g <- function(x) {
# return a vector of 1 and minus 1, since those are the values of the non-zero elements
return( c( rep( 1, 2*m ), rep( c(-1,1), m ) ) )
}
# ------------------------------------------------------------------
# rename lambda so it doesn't cause confusion with lambda in auxdata
eval_h <- function( x, obj_factor, hessian_lambda ) {
H <- t(A) %*% A
H <- unlist( lapply( 1:m, function(i) { H[i,1:i] } ) )
return( obj_factor * H )
}
eval_h_structure <- c( lapply( 1:m, function(x) { return( c(1:x) ) } ),
lapply( 1:m, function(x) { return( c() ) } ) )
# The starting point.
x0 = c( rep(0, m),
rep(1, m) )
# The constraint functions are bounded from below by zero.
constraint_lb = rep( 0, 2*m )
constraint_ub = rep( Inf, 2*m )
ipoptr_opts <- list( "jac_d_constant" = 'yes',
"hessian_constant" = 'yes',
"mu_strategy" = 'adaptive',
"max_iter" = 100,
"tol" = 1e-8 )
# Set up the auxiliary data.
auxdata <- new.env()
auxdata$m <- m
auxdata$A <- A
auxdata$y <- y
auxdata$lambda <- lambda
# COMPUTE SOLUTION WITH IPOPT.
# Compute the L1-regularized maximum likelihood estimator.
print( ipoptr( x0=x0,
eval_f=eval_f,
eval_grad_f=eval_grad_f,
eval_g=eval_g,
eval_jac_g=eval_jac_g,
eval_jac_g_structure=eval_jac_g_structure,
constraint_lb=constraint_lb,
constraint_ub=constraint_ub,
eval_h=eval_h,
eval_h_structure=eval_h_structure,
opts=ipoptr_opts,
ipoptr_environment=auxdata ) )
Mi punto es que si tienes nuevos datos, sólo necesitas
- actualización ( no reemplazar) la matriz de restricción y el vector de función objetivo para dar cuenta de las nuevas observaciones.
-
cambiar los puntos de partida del punto interior de
x0 = c( rep(0, m), rep(1, m) )
al vector de la solución que encontró anteriormente (antes de que se añadieran nuevos datos). La lógica aquí es la siguiente. Denota $ \beta_ {new}$ el nuevo vector de coeficientes (los que corresponden al conjunto de datos después de la actualización) y $ \beta_ {old}$ los originales. También denotan $ \beta_ {init}$ el vector x0 en el código de arriba (este es el inicio habitual para el método del punto interior). Entonces la idea es que si:
$$| \beta_ {init}- \beta_ {new}|_1>| \beta_ {new}- \beta_ {old}|_1 \quad (1)$$
entonces, uno puede conseguir $ \beta_ {new}$ mucho más rápido comenzando el punto interior desde $ \beta_ {old}$ en lugar de los ingenuos $ \beta_ {init}$ . La ganancia será tanto más importante cuando las dimensiones del conjunto de datos ( $n$ y $p$ ) son más grandes.
En cuanto a las condiciones en las que se mantiene la desigualdad (1), son:
- cuando $ \lambda $ es grande en comparación con $| \beta_ {OLS}|_1$ (este suele ser el caso cuando $p$ el número de variables de diseño es grande en comparación con $n$ el número de observaciones)
- cuando las nuevas observaciones no son patológicamente influyentes, por ejemplo cuando son coherentes con el proceso estocástico que ha generado los datos existentes.
- cuando el tamaño de la actualización es pequeño en relación con el tamaño de los datos existentes.
