Dlaczego Lars i Glmnet oferują różne rozwiązania problemu Lasso?

22

Chcę lepiej zrozumieć pakiety R Larsi Glmnetużywane do rozwiązania problemu Lasso: (dla zmiennych i próbek , patrz www.stanford.edu/~hastie/Papers/glmnet.pdf na stronie 3)

m ja n_{(β_{0} β) \in R^{p + 1}} [\frac{1}{2) N.} \sum_{ja = 1}^{N.} (y_{ja} - β_{0} - x_{ja}^{T.} β)^{2)} + λ | | β | |_{l_{1}}]

$min_{(\beta_0 \beta) \in R^{p+1}} \left[\frac{1}{2N}\sum_{i=1}^{N}(y_i-\beta_0-x_i^T\beta)^2 + \lambda||\beta ||_{l_{1}} \right]$

p

$p$

N

$N$

Dlatego zastosowałem je oba w tym samym zestawie danych zabawek. Niestety obie metody nie dają takich samych rozwiązań dla tych samych danych wejściowych. Czy ktoś ma pomysł, skąd ta różnica?

Otrzymałem wyniki w następujący sposób: Po wygenerowaniu niektórych danych (8 próbek, 12 cech, projekt Toeplitz, wszystko wyśrodkowane), obliczyłem całą ścieżkę Lasso za pomocą Larsa. Następnie uruchomiłem Glmnet, używając sekwencji lambd obliczonej przez Larsa (pomnożonej przez 0,5) i liczyłem na uzyskanie tego samego rozwiązania, ale nie zrobiłem tego.

Widać, że rozwiązania są podobne. Ale jak mogę wyjaśnić różnice? Mój kod znajduje się poniżej. Tutaj jest powiązane pytanie: GLMNET lub LARS do obliczania rozwiązań LASSO? , ale nie zawiera odpowiedzi na moje pytanie.

Ustawiać:

# Load packages.
library(lars)
library(glmnet)
library(MASS)

# Set parameters.
nb.features <- 12
nb.samples <- 8
nb.relevant.indices <- 3
snr <- 1
nb.lambdas <- 10

# Create data, not really important. 
sigma <- matrix(0, nb.features, nb.features)
for (i in (1:nb.features)) {
  for (j in (1:nb.features)) {
    sigma[i, j] <- 0.99 ^ (abs(i - j))
  }
}

x <- mvrnorm(n=nb.samples, rep(0, nb.features), sigma, tol=1e-6, empirical=FALSE)
relevant.indices <- sample(1:nb.features, nb.relevant.indices, replace=FALSE)
x <- scale(x)
beta <- rep(0, times=nb.features)
beta[relevant.indices] <- runif(nb.relevant.indices, 0, 1)
epsilon <- matrix(rnorm(nb.samples),nb.samples, 1)
simulated.snr <-(norm(x %*% beta, type="F")) / (norm(epsilon, type="F"))
epsilon <- epsilon * (simulated.snr / snr)
y <- x %*% beta + epsilon
y <- scale(y)

lars:

la <- lars(x, y, intercept=TRUE, max.steps=1000, use.Gram=FALSE)
co.lars <- as.matrix(coef(la, mode="lambda"))
print(round(co.lars, 4))

#          [,1] [,2] [,3]   [,4]   [,5]   [,6]    [,7]   [,8]    [,9]   [,10]
#  [1,]  0.0000    0    0 0.0000 0.0000 0.0000  0.0000 0.0000  0.0000  0.0000
#  [2,]  0.0000    0    0 0.0000 0.0000 0.1735  0.0000 0.0000  0.0000  0.0000
#  [3,]  0.0000    0    0 0.2503 0.0000 0.4238  0.0000 0.0000  0.0000  0.0000
#  [4,]  0.0000    0    0 0.1383 0.0000 0.7578  0.0000 0.0000  0.0000  0.0000
#  [5,] -0.1175    0    0 0.2532 0.0000 0.8506  0.0000 0.0000  0.0000  0.0000
#  [6,] -0.3502    0    0 0.2676 0.3068 0.9935  0.0000 0.0000  0.0000  0.0000
#  [7,] -0.4579    0    0 0.6270 0.0000 0.9436  0.0000 0.0000  0.0000  0.0000
#  [8,] -0.7848    0    0 0.9970 0.0000 0.9856  0.0000 0.0000  0.0000  0.0000
#  [9,] -0.3175    0    0 0.0000 0.0000 3.4488  0.0000 0.0000 -2.1714  0.0000
# [10,] -0.4842    0    0 0.0000 0.0000 4.7731  0.0000 0.0000 -3.4102  0.0000
# [11,] -0.4685    0    0 0.0000 0.0000 4.7958  0.0000 0.1191 -3.6243  0.0000
# [12,] -0.4364    0    0 0.0000 0.0000 5.0424  0.0000 0.3007 -4.0694 -0.4903
# [13,] -0.4373    0    0 0.0000 0.0000 5.0535  0.0000 0.3213 -4.1012 -0.4996
# [14,] -0.4525    0    0 0.0000 0.0000 5.6876 -1.5467 1.5095 -4.7207  0.0000
# [15,] -0.4593    0    0 0.0000 0.0000 5.7355 -1.6242 1.5684 -4.7440  0.0000
# [16,] -0.4490    0    0 0.0000 0.0000 5.8601 -1.8485 1.7767 -4.9291  0.0000
#         [,11]  [,12]
#  [1,]  0.0000 0.0000
#  [2,]  0.0000 0.0000
#  [3,]  0.0000 0.0000
#  [4,] -0.2279 0.0000
#  [5,] -0.3266 0.0000
#  [6,] -0.5791 0.0000
#  [7,] -0.6724 0.2001
#  [8,] -1.0207 0.4462
#  [9,] -0.4912 0.1635
# [10,] -0.5562 0.2958
# [11,] -0.5267 0.3274
# [12,]  0.0000 0.2858
# [13,]  0.0000 0.2964
# [14,]  0.0000 0.1570
# [15,]  0.0000 0.1571

glmnet z lambda = (lambda_lars / 2):

glm2 <- glmnet(x, y, family="gaussian", lambda=(0.5 * la$lambda), thresh=1e-16)
co.glm2 <- as.matrix(t(coef(glm2, mode="lambda")))
print(round(co.glm2, 4))

#     (Intercept)      V1 V2 V3     V4     V5     V6      V7     V8      V9
# s0            0  0.0000  0  0 0.0000 0.0000 0.0000  0.0000 0.0000  0.0000
# s1            0  0.0000  0  0 0.0000 0.0000 0.0000  0.0000 0.0000  0.0000
# s2            0  0.0000  0  0 0.2385 0.0000 0.4120  0.0000 0.0000  0.0000
# s3            0  0.0000  0  0 0.2441 0.0000 0.4176  0.0000 0.0000  0.0000
# s4            0  0.0000  0  0 0.2466 0.0000 0.4200  0.0000 0.0000  0.0000
# s5            0  0.0000  0  0 0.2275 0.0000 0.4919  0.0000 0.0000  0.0000
# s6            0  0.0000  0  0 0.1868 0.0000 0.6132  0.0000 0.0000  0.0000
# s7            0 -0.2651  0  0 0.2623 0.1946 0.9413  0.0000 0.0000  0.0000
# s8            0 -0.6609  0  0 0.7328 0.0000 1.6384  0.0000 0.0000 -0.5755
# s9            0 -0.4633  0  0 0.0000 0.0000 4.6069  0.0000 0.0000 -3.2547
# s10           0 -0.4819  0  0 0.0000 0.0000 4.7546  0.0000 0.0000 -3.3929
# s11           0 -0.4767  0  0 0.0000 0.0000 4.7839  0.0000 0.0567 -3.5122
# s12           0 -0.4715  0  0 0.0000 0.0000 4.7915  0.0000 0.0965 -3.5836
# s13           0 -0.4510  0  0 0.0000 0.0000 5.6237 -1.3909 1.3898 -4.6583
# s14           0 -0.4552  0  0 0.0000 0.0000 5.7064 -1.5771 1.5326 -4.7298
#         V10     V11    V12
# s0   0.0000  0.0000 0.0000
# s1   0.0000  0.0000 0.0000
# s2   0.0000  0.0000 0.0000
# s3   0.0000  0.0000 0.0000
# s4   0.0000  0.0000 0.0000
# s5   0.0000 -0.0464 0.0000
# s6   0.0000 -0.1293 0.0000
# s7   0.0000 -0.4868 0.0000
# s8   0.0000 -0.8803 0.3712
# s9   0.0000 -0.5481 0.2792
# s10  0.0000 -0.5553 0.2939
# s11  0.0000 -0.5422 0.3108
# s12  0.0000 -0.5323 0.3214
# s13 -0.0503  0.0000 0.1711
# s14  0.0000  0.0000 0.1571

— Andre
źródło

20

Wreszcie udało nam się stworzyć to samo rozwiązanie za pomocą obu metod! Pierwszy problem polega na tym, że glmnet rozwiązuje problem lasso, jak stwierdzono w pytaniu, ale lars ma nieco inną normalizację w funkcji celu, zastępuje przez . Po drugie, obie metody normalizują dane w różny sposób, więc normalizację należy wyłączyć podczas wywoływania metod. $\frac{1}{2N}$ $\frac{1}{2}$

Aby to odtworzyć i przekonać się, że te same rozwiązania problemu lasso można obliczyć za pomocą lars i glmnet, należy zmienić następujące wiersze w powyższym kodzie:

la <- lars(X,Y,intercept=TRUE, max.steps=1000, use.Gram=FALSE)

do

la <- lars(X,Y,intercept=TRUE, normalize=FALSE, max.steps=1000, use.Gram=FALSE)

i

glm2 <- glmnet(X,Y,family="gaussian",lambda=0.5*la$lambda,thresh=1e-16)

do

glm2 <- glmnet(X,Y,family="gaussian",lambda=1/nbSamples*la$lambda,standardize=FALSE,thresh=1e-16)

— Andre
źródło

1

Cieszę się, że to rozgryzłeś. Czy są jakieś przemyślenia, która metoda normalizacji ma większy sens? Właściwie uzyskałem gorsze wyniki dzięki normalizacji w glmnet (dla lasso) i nadal nie jestem pewien dlaczego.

— Ben Ogorek

Właściwie normalizuję dane od razu i stosuję te metody i porównuję je, jeśli są podobne. Zmienne o mniejszych efektach zwykle mają różne współczynniki

— KarthikS

0

Oczywiście, jeśli metody wykorzystują różne modele, otrzymasz różne odpowiedzi. Odjęcie warunków przechwytywania nie prowadzi do modelu bez przechwytywania, ponieważ zmienią się współczynniki najlepszego dopasowania i nie zmienisz ich w miarę zbliżania się do niego. Musisz dopasować ten sam model do obu metod, jeśli chcesz uzyskać te same lub prawie takie same odpowiedzi.

— Michael R. Chernick
źródło

1

Tak, masz rację, metody używają nieco innych modeli, nie byłem tego świadomy. Dzięki za podpowiedź. (Różnice wyjaśnię bardziej szczegółowo w osobnej odpowiedzi)

— Andre

-2

Wyniki muszą być takie same. pakiet lars używa domyślnie type = "lar", zmień tę wartość na type = "lasso". Po prostu obniż parametr „thresh = 1e-16” dla glmnet, ponieważ opadanie współrzędnych opiera się na zbieżności.

— Marcool Lopez Cruz
źródło

2

Dziękuję za Twoją odpowiedź. Być może źle to interpretuję, ale wydaje się to sprzeczne z rezolucją opublikowaną w odpowiedzi Andre sześć lat temu. Rozważ rozwinięcie swojego postu, tak aby zawierało pełniejsze wyjaśnienie tego, co próbujesz powiedzieć, i wskazanie, dlaczego uważamy, że jest ono poprawne, a drugie nie.

— whuber