套索是通过 LARS 拟合的（一个迭代过程，从一些初始估计开始）。默认情况下但您可以在大多数实现中更改它（并用您已经拥有最接近是，你将不得不步入的 LARS 迭代次数越少。$\beta^0$$\beta^0=0_p$$\beta^*_{old}$$\beta^*_{old}$$\beta_{new}^*$$\beta_{new}^*$

编辑：

由于来自 user2763361我的评论，我在我的原始答案中添加了更多细节。

从下面的评论中，我收集到 user2763361 建议补充我的原始答案，将其变成可以直接使用（下架）同时也非常有效的答案。

为了完成第一部分，我将通过一个玩具示例逐步说明我提出的解决方案。为了满足第二部分，我将使用最近的高质量内点求解器来完成。这是因为，使用可以通过内点方法解决套索问题的库，而不是尝试破解 LARS 或单纯形算法以从非标准起点（尽管第二个地点也是可能的）。

请注意，有时（在旧书中）声称求解线性规划的内点法比单纯形法慢，这在很久以前可能是正确的，但今天通常不是正确的，对于大规模问题当然也不是正确的（这就是为什么大多数专业图书馆都喜欢cplex使用内点算法）并且问题至少隐含地涉及大规模问题。另请注意，我使用的内点求解器完全处理稀疏矩阵，因此我认为与 LARS 不会有很大的性能差距（使用 LARS 的最初动机是当时许多流行的 LP 求解器不能很好地处理稀疏矩阵，并且这些是 LASSO 问题的特征）。

ipopt库中的一个（非常）好的内点算法开源实现COIN-OR。我将使用的另一个原因ipopt是它有一个 R 接口，ipoptr. 你会在这里找到更详尽的安装指南，下面我给出了安装它的标准命令ubuntu。

在中bash，做：

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

然后，以 root 身份R执行（我假设svn已将 subversion 文件复制到~/默认情况下）：

install.packages("~/CoinIpopt/Ipopt/contrib/RInterface",repos=NULL,type="source")

从这里，我给出一个小例子（主要来自 Jelmer Ypma 作为他的包装器的一部分给出的玩具示例 R）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 ) )

我的意思是，如果你有新数据，你只需要

1. 更新（而不是替换）约束矩阵和目标函数向量以解释新的观察结果。

2. 将内部点的起点从

   x0 = c（代表（0，米），代表（1，米））

   到您之前找到的解决方案向量（在添加新数据之前）。这里的逻辑如下。表示新的系数向量（对应于更新后的数据集）和原始向量。还要在上面的代码中表示向量（这是内点方法的通常开始）。那么这个想法是，如果：$\beta_{new}$$\beta_{old}$$\beta_{init}$x0

然后，而不是幼稚的开始内点，可以更快 } 。和）的维度更大时，增益将变得更加重要。$\beta_{new}$$\beta_{old}$$\beta_{init}$$n$$p$

至于不等式（1）成立的条件，分别是：

• 当与相比较大时（这通常是当与观察数量相比较大时）$\lambda$$|\beta_{OLS}|_1$$p$$n$
• 当新的观察结果没有病理影响时，例如当它们与生成现有数据的随机过程一致时。
• 当更新的大小相对于现有数据的大小较小时。