Replizieren der Ergebnisse für die lineare glmnet-Regression mithilfe eines generischen Optimierers

Wie der Titel schon sagt, versuche ich, die Ergebnisse von glmnet linear mit dem LBFGS-Optimierer aus der Bibliothek zu replizieren lbfgs. Mit diesem Optimierer können wir einen L1-Regularisierungsbegriff hinzufügen, ohne uns um die Differenzierbarkeit kümmern zu müssen, solange unsere Zielfunktion (ohne den L1-Regularisierungsbegriff) konvex ist.

Das Problem der linearen Regression des elastischen Netzes im glmnet-Papier ist gegeben durch wobei ist die Entwurfsmatrix, ist der Beobachtungsvektor, ist der elastische Netzparameter und ist der Regularisierungsparameter. Der Operator bezeichnet die übliche Lp-Norm.

min_{β \in R^{p}} \frac{1}{2 n} ‖ β_{0} + X β - y ‖_{2}^{2} + α λ ‖ β ‖_{1} + \frac{1}{2} (1 - α) λ ‖ β ‖_{2}^{2}

$\min_{\beta \in \mathbb{R}^p} \frac{1}{2n}\Vert \beta_0 + X\beta - y \Vert_2^2 + \alpha \lambda \Vert \beta\Vert_1 + \frac{1}{2}(1-\alpha)\lambda\Vert\beta\Vert^2_2$

X \in R^{n \times p}

$X \in \mathbb{R}^{n \times p}$

y \in R^{p}

$y \in \mathbb{R}^p$

α \in [0, 1]

$\alpha \in [0,1]$

λ > 0

$\lambda > 0$

‖ x ‖_{p}

$\Vert x \Vert_p$

Der folgende Code definiert die Funktion und enthält dann einen Test zum Vergleichen der Ergebnisse. Wie Sie sehen können, sind die Ergebnisse akzeptabel, wenn alpha = 1, aber weit entfernt von den Werten von alpha < 1.Der Fehler wird schlimmer, wenn wir von alpha = 1bis gehen alpha = 0, wie das folgende Diagramm zeigt (die "Vergleichsmetrik" ist der mittlere euklidische Abstand zwischen den Parameterschätzungen von glmnet und lbfgs für einen gegebenen Regularisierungspfad).

Okay, hier ist der Code. Ich habe wo immer möglich Kommentare hinzugefügt. Meine Frage ist: Warum unterscheiden sich meine Ergebnisse von denen glmnetfür Werte von alpha < 1? Es hat eindeutig mit dem L2-Regularisierungsbegriff zu tun, aber soweit ich das beurteilen kann, habe ich diesen Begriff genau gemäß dem Papier implementiert. Jede Hilfe wäre sehr dankbar!

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 Ich habe Ihren Code ausprobiert. Hier sind einige Tests (ich habe einige kleinere Änderungen vorgenommen, um sie an die Struktur von glmnet anzupassen. Beachten Sie, dass wir den Intercept-Term nicht regulieren und die Verlustfunktionen skaliert werden müssen.) Dies ist für alpha = 0, aber Sie können jeden ausprobieren alpha- die Ergebnisse stimmen nicht überein.

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)))

r regression regularization glmnet elastic-net user3294195
quelle

Ich bin mir nicht sicher, ob Ihre Frage zum Thema gehört (ich denke, es könnte sein, dass es sich um die zugrunde liegende Optimierungstechnik handelt), und ich kann Ihren Code jetzt nicht wirklich überprüfen, aber ich lbfgsmöchte einen Punkt zum orthantwise_cArgument der glmnetÄquivalenz ansprechen.

Firebug

Das Problem ist nicht wirklich mit lbfgsund orthantwise_cwie damals alpha = 1ist die Lösung fast genau die gleiche mit glmnet. Es hat mit der L2-Regularisierungsseite der Dinge zu tun, dh wann alpha < 1. Ich denke, eine Änderung an der Definition von vorzunehmen SSEund diese SSE_grzu beheben, aber ich bin mir nicht sicher, wie die Änderung aussehen soll - soweit ich weiß, sind diese Funktionen genau wie im glmnet-Dokument beschrieben definiert.

user3294195

Dies kann eher eine Stapelüberlauf- oder Programmierfrage sein.

Matthew Gunn

Ich dachte, es hat mehr mit Optimierung und Regularisierung zu tun als mit dem Code selbst, weshalb ich ihn hier gepostet habe.

user3294195

Für eine reine Optimierungsfrage ist scicomp.stackexchange.com ebenfalls eine Option. Ich bin mir nicht sicher, ob dort sprachspezifische Fragen zum Thema stehen. (zB "mach das in R")

GeoMatt22

tl; dr version:

Das Ziel enthält implizit einen Skalierungsfaktor , wobei die Standardabweichung der Stichprobe ist. $\hat{s} = sd(y)$ $sd(y)$

Längere Version

Wenn Sie das Kleingedruckte der glmnet-Dokumentation lesen, sehen Sie:

Beachten Sie, dass die Zielfunktion für "Gauß" ist
               1/2  RSS/nobs + lambda*penalty,                  
und für die anderen Modelle ist es
               -loglik/nobs + lambda*penalty.                   
Beachten Sie auch, dass für '"Gauß'sches"' glmnet 'y standardisiert, um eine Einheitsvarianz zu haben, bevor seine Lambda-Sequenz berechnet wird (und dann die resultierenden Koeffizienten nicht standardisiert); Wenn Sie Ergebnisse mit anderer Software reproduzieren / vergleichen möchten, geben Sie am besten ein standardisiertes y an.

Dies bedeutet nun, dass das Ziel tatsächlich und dieses glmnet meldet .

\frac{1}{2 n} ‖ y / \hat{s} - X β ‖_{2}^{2} + λ α ‖ β ‖_{1} + λ (1 - α) ‖ β ‖_{2}^{2},

$\frac{1}{2n} \lVert y/\hat{s}-X\beta \rVert_2^2 + \lambda\alpha \lVert \beta \rVert_1 + \lambda(1-\alpha)\lVert \beta \rVert_2^2,$

\tilde{β} = \hat{s} β

$\tilde{\beta} = \hat{s} \beta$

Wenn Sie nun ein reines Lasso ( ) verwenden, bedeutet eine Nichtstandardisierung von glmnets , dass die Antworten gleichwertig sind. Auf der anderen Seite müssen Sie bei einem reinen Grat die Strafe um einen Faktor skalieren , damit der Pfad übereinstimmt, da ein zusätzlicher Faktor von aus dem Quadrat herausspringt in der Strafe. Für Intermediate gibt es keine einfache Möglichkeit, die Strafe für Koeffizienten zu skalieren, um die Ausgabe zu reproduzieren . $\alpha=1$ $\tilde{\beta}$ $1/\hat{s}$ glmnet $\hat{s}$ $\ell_2$ $\alpha$ glmnets

Sobald ich das skaliere, um eine Einheitsvarianz zu erhalten, finde ich $y$

was immer noch nicht genau passt. Dies scheint auf zwei Dinge zurückzuführen zu sein:

Die Lambda-Sequenz ist möglicherweise zu kurz, als dass der zyklische Koordinatenabstiegsalgorithmus für den Warmstart vollständig konvergent wäre.
Ihre Daten enthalten keinen Fehlerterm ( der Regression ist 1). $R^2$
Beachten Sie auch, dass der bereitgestellte Code einen Fehler enthält, der lambda[2]für die anfängliche Anpassung erforderlich ist. Dies sollte jedoch der Fall sein lambda[1].

Sobald ich die Punkte 1-3 korrigiere, erhalte ich das folgende Ergebnis (obwohl YMMV vom zufälligen Startwert abhängt):

Andrew M.
quelle

Replizieren der Ergebnisse für die lineare glmnet-Regression mithilfe eines generischen Optimierers

Antworten:

tl; dr version:

Längere Version