# Updating the lasso fit with new observations

I am fitting an L1-regularized linear regression to a very large dataset (with n>>p.) The variables are known in advance, but the observations arrive in small chunks. I would like to maintain the lasso fit after each chunk.

I can obviously re-fit the entire model after seeing each new set of observations. This, however, would be pretty inefficient given that there is a lot of data. The amount of new data that arrives at each step is very small, and the fit is unlikely to change much between steps.

Is there anything I can do to reduce the overall computational burden?

I was looking at the LARS algorithm of Efron et al., but would be happy to consider any other fitting method if it can be made to "warm-start" in the way described above.

Notes:

1. I am mainly looking for an algorithm, but pointers to existing software packages that can do this may also prove insightful.
2. In addition to the current lasso trajectories, the algorithm is of course welcome to keep other state.

Bradley Efron, Trevor Hastie, Iain Johnstone and Robert Tibshirani, Least Angle Regression, Annals of Statistics (with discussion) (2004) 32(2), 407--499.

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) )
}
# ------------------------------------------------------------------
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',
"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_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

1. update (not replace) the constraint matrix and objective function vector to account for the new observations.
2. 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.
• Thanks, but I am afraid I don't follow. LARS produces a piecewise-linear path (with exactly $p+1$ points for the least angles and possibly more points for the lasso.) Each point has its own set of $\beta$. When we add more observations, all the betas can move (except $\beta^0$, which is always $0_p$.) Please could you expand on your answer? Thanks. – NPE Nov 17 '10 at 11:34
• @aix: do you want to update the whole lasso path or just the solution ? (i.e. is the sparsity penalty fixed?). – user603 Nov 17 '10 at 16:57
• I was looking to update the entire path. However, if there's a good way to do it for a fixed penalty ($\lambda$ in the formula below), this may be a good start. Is this what you are proposing? $$\hat{\beta}^{lasso} = \underset{\beta}{\operatorname{argmin}} \left \{ {1 \over 2} \sum_{i=1}^N(y_i-\beta_0-\sum_{j=1}^p x_{ij} \beta_j)^2 + \lambda \sum_{j=1}^p |\beta_j| \right \}$$ – NPE Nov 17 '10 at 17:48
• @aix. Yes, it all depends on the implementation you use and the facilities you have access to. For example: if you have access to a good lp solver, you can feed it with the past optimal values of $\beta$ and it'll carry the 1-2 step to the new solution very efficiently. You should add these details to your question. – user603 Nov 20 '10 at 16:27
• Any libraries you know that can do this without needing to edit the source code? – user2763361 Jun 6 '14 at 3:19