Replicando resultados para regressão linear glmnet usando um otimizador genérico

Como o título indica, estou tentando replicar os resultados do glmnet linear usando o otimizador LBFGS da biblioteca lbfgs. Esse otimizador nos permite adicionar um termo regularizador L1 sem ter que nos preocupar com diferenciabilidade, desde que nossa função objetivo (sem o termo regularizador L1) seja convexa.

O problema da regressão linear líquida elástica no papel glmnet é dado por que X \ in \ mathbb {R} ^ {n \ times p} é a matriz de design, y \ in \ mathbb {R} ^ p é o vetor de observações, \ alpha \ in [0,1] é o parâmetro da rede elástica e \ lambda> 0 é o parâmetro de regularização. O operador \ Vert x \ Vert_p denota a norma Lp usual.

O código abaixo define a função e inclui um teste para comparar os resultados. Como você pode ver, os resultados são aceitáveis ​​quando alpha = 1, mas estão longe dos valores de alpha < 1.O erro piora à medida que avançamos alpha = 1para alpha = 0, como mostra o gráfico a seguir (a "métrica de comparação" é a distância euclidiana média entre as estimativas de parâmetros do glmnet e lbfgs para um determinado caminho de regularização).

Ok, então aqui está o código. Adicionei comentários sempre que possível. Minha pergunta é: por que meus resultados são diferentes dos glmnetvalores de alpha < 1? Claramente, isso tem a ver com o termo de regularização L2, mas, tanto quanto posso dizer, implementei esse termo exatamente conforme o documento. Qualquer ajuda seria muito apreciada!

library(lbfgs)
linreg_lbfgs <- function(X, y, alpha = 1, scale = TRUE, lambda) {
  p <- ncol(X) + 1; n <- nrow(X); nlambda <- length(lambda)

  # Scale design matrix
  if (scale) {
    means <- colMeans(X)
    sds <- apply(X, 2, sd)
    sX <- (X - tcrossprod(rep(1,n), means) ) / tcrossprod(rep(1,n), sds)
  } else {
    means <- rep(0,p-1)
    sds <- rep(1,p-1)
    sX <- X
  }
  X_ <- cbind(1, sX)

  # loss function for ridge regression (Sum of squared errors plus l2 penalty)
  SSE <- function(Beta, X, y, lambda0, alpha) {
    1/2 * (sum((X%*%Beta - y)^2) / length(y)) +
      1/2 * (1 - alpha) * lambda0 * sum(Beta[2:length(Beta)]^2) 
                    # l2 regularization (note intercept is excluded)
  }

  # loss function gradient
  SSE_gr <- function(Beta, X, y, lambda0, alpha) {
    colSums(tcrossprod(X%*%Beta - y, rep(1,ncol(X))) *X) / length(y) + # SSE grad
  (1-alpha) * lambda0 * c(0, Beta[2:length(Beta)]) # l2 reg grad
  }

  # matrix of parameters
  Betamat_scaled <- matrix(nrow=p, ncol = nlambda)

  # initial value for Beta
  Beta_init <- c(mean(y), rep(0,p-1)) 

  # parameter estimate for max lambda
  Betamat_scaled[,1] <- lbfgs(call_eval = SSE, call_grad = SSE_gr, vars = Beta_init, 
                              X = X_, y = y, lambda0 = lambda[2], alpha = alpha,
                              orthantwise_c = alpha*lambda[2], orthantwise_start = 1, 
                              invisible = TRUE)$par

  # parameter estimates for rest of lambdas (using warm starts)
  if (nlambda > 1) {
    for (j in 2:nlambda) {
      Betamat_scaled[,j] <- lbfgs(call_eval = SSE, call_grad = SSE_gr, vars = Betamat_scaled[,j-1], 
                                  X = X_, y = y, lambda0 = lambda[j], alpha = alpha,
                                  orthantwise_c = alpha*lambda[j], orthantwise_start = 1, 
                                  invisible = TRUE)$par
    }
  }

  # rescale Betas if required
  if (scale) {
    Betamat <- rbind(Betamat_scaled[1,] -
colSums(Betamat_scaled[-1,]*tcrossprod(means, rep(1,nlambda)) / tcrossprod(sds, rep(1,nlambda)) ), Betamat_scaled[-1,] / tcrossprod(sds, rep(1,nlambda)) )
  } else {
    Betamat <- Betamat_scaled
  }
  colnames(Betamat) <- lambda
  return (Betamat)
}

# CODE FOR TESTING
# simulate some linear regression data
n <- 100
p <- 5
X <- matrix(rnorm(n*p),n,p)
true_Beta <- sample(seq(0,9),p+1,replace = TRUE)
y <- drop(cbind(1,X) %*% true_Beta)

library(glmnet)

# function to compare glmnet vs lbfgs for a given alpha
glmnet_compare <- function(X, y, alpha) {
  m_glmnet <- glmnet(X, y, nlambda = 5, lambda.min.ratio = 1e-4, alpha = alpha)
  Beta1 <- coef(m_glmnet)
  Beta2 <- linreg_lbfgs(X, y, alpha = alpha, scale = TRUE, lambda = m_glmnet$lambda)
  # mean Euclidean distance between glmnet and lbfgs results
  mean(apply (Beta1 - Beta2, 2, function(x) sqrt(sum(x^2))) ) 
}

# compare results
alpha_seq <- seq(0,1,0.2)
plot(alpha_seq, sapply(alpha_seq, function(alpha) glmnet_compare(X,y,alpha)), type = "l", ylab = "Comparison metric")

@ hxd1011 Tentei o seu código, eis alguns testes (fiz alguns pequenos ajustes para combinar com a estrutura do glmnet - observe que não regularizamos o termo de interceptação e as funções de perda devem ser dimensionadas). Isto é para alpha = 0, mas você pode tentar qualquer um alpha- os resultados não correspondem.

rm(list=ls())
set.seed(0)
# simulate some linear regression data
n <- 1e3
p <- 20
x <- matrix(rnorm(n*p),n,p)
true_Beta <- sample(seq(0,9),p+1,replace = TRUE)
y <- drop(cbind(1,x) %*% true_Beta)

library(glmnet)
alpha = 0

m_glmnet = glmnet(x, y, alpha = alpha, nlambda = 5)

# linear regression loss and gradient
lr_loss<-function(w,lambda1,lambda2){
  e=cbind(1,x) %*% w -y
  v= 1/(2*n) * (t(e) %*% e) + lambda1 * sum(abs(w[2:(p+1)])) + lambda2/2 * crossprod(w[2:(p+1)])
  return(as.numeric(v))
}

lr_loss_gr<-function(w,lambda1,lambda2){
  e=cbind(1,x) %*% w -y
  v= 1/n * (t(cbind(1,x)) %*% e) + c(0, lambda1*sign(w[2:(p+1)]) + lambda2*w[2:(p+1)])
  return(as.numeric(v))
}

outmat <- do.call(cbind, lapply(m_glmnet$lambda, function(lambda) 
  optim(rnorm(p+1),lr_loss,lr_loss_gr,lambda1=alpha*lambda,lambda2=(1-alpha)*lambda,method="L-BFGS")$par
))

glmnet_coef <- coef(m_glmnet)
apply(outmat - glmnet_coef, 2, function(x) sqrt(sum(x^2)))

tl; dr versão:

O objetivo contém implicitamente um fator de escala , em que é o desvio padrão da amostra.$\hat{s} = sd(y)$$sd(y)$

Versão mais longa

Se você ler a cópia fina da documentação do glmnet, verá:

Observe que a função objetivo para '"gaussian"' é

               1/2  RSS/nobs + lambda*penalty,                  

e para os outros modelos é

               -loglik/nobs + lambda*penalty.                   

Observe também que para '"gaussian"', 'glmnet' padroniza y para ter variação unitária antes de calcular sua sequência lambda (e depois padronizar os coeficientes resultantes); se você deseja reproduzir / comparar resultados com outro software, é melhor fornecer um y padronizado.

Agora, isso significa que o objetivo é, na verdade, e que glmnet reporta .

Agora, quando você estava usando um laço puro ( ), a não padronização do da glmnet significa que as respostas são equivalentes. Por outro lado, com uma cordilheira pura, é necessário escalar a penalidade por um fator para que o caminho seja aceito, porque um fator extra de sai da praça na penalidade . Para intermediário , não há uma maneira fácil de escalar a penalidade de coeficientes para reproduzir a saída.$\alpha=1$$\tilde{\beta}$$1/\hat{s}$glmnet$\hat{s}$$\ell_2$$\alpha$glmnets

Depois de escalar para ter variação de unidade, encontro $y$

que ainda não corresponde exatamente. Isso parece dever-se a duas coisas:

1. A sequência lambda pode ser muito curta para que o algoritmo de descida de coordenadas cíclicas de inicialização a quente seja totalmente convergente.
2. Não há termo de erro nos seus dados (o da regressão é 1).$R^2$
3. Observe também que existe um erro no código, conforme o necessário lambda[2]para o ajuste inicial, mas deve ser lambda[1].

Depois de corrigir os itens 1-3, recebo o seguinte resultado (embora o YMMV dependa da semente aleatória):