The lasso is fitted through LARS (an iterative process, that starts at some initial estimate $\beta^0$). By default $\beta^0=0_p$ but you can change this in most implementation (and replace it by the optimal $\beta^*_{old}$ you already have). The closest $\beta^*_{old}$ is to $\beta_{new}^*$, the smaller the number of LARS iteration you will have to step to get to $\beta_{new}^*$.
EDIT:
Due to the comments from user2763361 I add more details to my original answer.
From the comments below I gather that user2763361 suggests to complement my original answer to turn it into one that can be used directly (off the shelves) while also being very efficient.
To do the first part, I will illustrate the solution I propose step by step on a toy example. To satisfy the second part, I will do so using a recent, high quality interior point solver.
This is because, it is easier to obtain an high performance implementation of the solution I propose using a library that can solve the lasso problem by the interior point approach rather than trying to hack the LARS or simplex algorithm to start the optimization from a non-standard starting point (though that second venue is also possible).
Note that it is sometimes claimed (in older books) that the interior point approach to solving linear programs is slower than the simplex approach and that may have been true a long time ago but it's generally not true today and certainly not true for large scale problems (this is why most professional libraries like cplex use the interior point algorithm) and the question is at least implicitly about large scale problems.
Note also that the interior point solver I use fully handles sparse matrices so I don t think there will be a large performance gap with LARS (an original motivation for using LARS was that many popular LP solvers at the time were not handling sparse matrices well and these are a characteristic features of the LASSO problem).
A (very) good open source implementation of the interior point algorithm is ipopt, in the COIN-OR library. Another reason I will be using ipopt is that it has has an R interface, ipoptr. You will find more exhaustive installation guide here, below I give the standard commands to install it in ubuntu.
in the bash, do:
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
Then, as root, in R do (I assume svn has copied the subversion file in ~/ as it does by default):
install.packages("~/CoinIpopt/Ipopt/contrib/RInterface",repos=NULL,type="source")
From here, I'm giving a small example (mostly from the toy example given by Jelmer Ypma as part of his R wraper to 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 ) )
My point is, if you have new data in, you just need to
- update (not replace) the constraint matrix and objective function vector to account for the new observations.
change the starting points of the interior point from
x0 = c( rep(0, m),
rep(1, m) )
to the vector of solution you found previously (before new data was added in). The logic here is as follows. Denote $\beta_{new}$ the new vector of coefficients (the ones corresponding to the data set after the update) and $\beta_{old}$ the original ones. Also denote $\beta_{init}$ the vector x0 in the code above (this is the usual start for the interior point method). Then the idea is that if:
$$|\beta_{init}-\beta_{new}|_1>|\beta_{new}-\beta_{old}|_1\quad(1)$$
then, one can get $\beta_{new}$ much faster by starting the interior point from
$\beta_{old}$ rather than the naive $\beta_{init}$. The gain will be all the more important when the dimensions of the data set ($n$ and $p$) are larger.
As for the conditions under which inequality (1) holds, they are:
- when $\lambda$ is large compared to $|\beta_{OLS}|_1$ (this is usually the case when $p$, the number of design variables is large compared to $n$, the number of observations)
- when the new observations are not pathologically influential, e.g. for example when they are consistent with the stochastic process that has generated the existing data.
- when the size of the update is small relative to the size of the existing data.
