Warum konvergiert Lasso nicht bei einem Bestrafungsparameter?

7

Um herauszufinden, wie die LASSORegression funktioniert, habe ich einen kleinen Code geschrieben, der die LASSORegression durch Auswahl des besten Alpha-Parameters optimieren soll .

Ich kann nicht herausfinden, warum die LASSORegression nach der Kreuzvalidierung so instabile Ergebnisse für den Alpha-Parameter liefert.

Hier ist mein Python-Code:

from sklearn.linear_model import Lasso
from sklearn.cross_validation import KFold
from matplotlib import pyplot as plt

# generate some sparse data to play with
import numpy as np
import pandas as pd 
from scipy.stats import norm
from scipy.stats import uniform

### generate your own data here

n = 1000

x1x2corr = 1.1
x1x3corr = 1.0
x1 = range(n) + norm.rvs(0, 1, n) + 50
x2 =  map(lambda aval: aval*x1x2corr, x1) + norm.rvs(0, 2, n) + 500
y = x1 + x2 #+ norm.rvs(0,10, n)

Xdf = pd.DataFrame()
Xdf['x1'] = x1
Xdf['x2'] = x2

X = Xdf.as_matrix()

# Split data in train set and test set
n_samples = X.shape[0]
X_train, y_train = X[:n_samples / 2], y[:n_samples / 2]
X_test, y_test = X[n_samples / 2:], y[n_samples / 2:]

kf = KFold(X_train.shape[0], n_folds = 10, )
alphas = np.logspace(-16, 8, num = 1000, base = 2)

e_alphas = list()
e_alphas_r = list()  # holds average r2 error
for alpha in alphas:
    lasso = Lasso(alpha=alpha, tol=0.004)
    err = list()
    err_2 = list()
    for tr_idx, tt_idx in kf:
        X_tr, X_tt = X_train[tr_idx], X_test[tt_idx]
        y_tr, y_tt = y_train[tr_idx], y_test[tt_idx]
        lasso.fit(X_tr, y_tr)
        y_hat = lasso.predict(X_tt)

        # returns the coefficient of determination (R^2 value)
        err_2.append(lasso.score(X_tt, y_tt))

        # returns MSE
        err.append(np.average((y_hat - y_tt)**2))
    e_alphas.append(np.average(err))
    e_alphas_r.append(np.average(err_2))

## print out the alpha that gives the minimum error
print 'the minimum value of error is ', e_alphas[e_alphas.index(min(e_alphas))]
print ' the minimizer is ',  alphas[e_alphas.index(min(e_alphas))]

##  <<< plotting alphas against error >>>

plt.figsize = (15, 15)
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(alphas, e_alphas, 'b-')
ax.plot(alphas, e_alphas_r, 'g--')
ax.set_ylim(min(e_alphas),max(e_alphas))
ax.set_xlim(min(alphas),max(alphas))
ax.set_xlabel("alpha")
plt.show()

Wenn Sie diesen Code wiederholt ausführen, ergeben sich für Alpha völlig unterschiedliche Ergebnisse:

>>> 
the minimum value of error is  3.99254192539
 the minimizer is  1.52587890625e-05
>>> ================================ RESTART ================================
>>> 
the minimum value of error is  4.07412455842
 the minimizer is  6.45622425334
>>> ================================ RESTART ================================
>>> 
the minimum value of error is  4.25898253597
 the minimizer is  1.52587890625e-05
>>> ================================ RESTART ================================
>>> 
the minimum value of error is  3.79392968781
 the minimizer is  28.8971008254
>>> 

Warum konvergiert der Alpha-Wert nicht richtig? Ich weiß, dass meine Daten synthetisch sind, aber die Verteilung ist dieselbe. Auch die Variation ist in x1und sehr gering x2.

Was könnte dazu führen, dass dies so instabil ist?

Dasselbe, was in R geschrieben ist, führt zu unterschiedlichen Ergebnissen - es gibt immer den höchstmöglichen Wert für Alpha als "optimal_alpha" zurück.

Ich habe das auch in R geschrieben, was mir eine etwas andere Antwort gibt, die ich nicht weiß warum?

library(glmnet)
library(lars)
library(pracma)

set.seed(1)
k = 2 # number of features selected 

n = 1000

x1x2corr = 1.1
x1 = seq(n) + rnorm(n, 0, 1) + 50
x2 =  x1*x1x2corr + rnorm(n, 0, 2) + 500
y = x1 + x2 

filter_out_label <- function(col) {col!="y"}

alphas = logspace(-5, 6, 100)

for (alpha in alphas){
  k = 10
  optimal_alpha = NULL
  folds <- cut(seq(1, nrow(df)), breaks=k, labels=FALSE)
  total_mse = 0
  min_mse = 10000000
  for(i in 1:k){
    # Segement your data by fold using the which() function
    testIndexes <- which(folds==i, arr.ind=TRUE)
    testData <- df[testIndexes, ]
    trainData <- df[-testIndexes, ]

    fit <- lars(as.matrix(trainData[Filter(filter_out_label, names(df))]),
                trainData$y,
                type="lasso")
    # predict
    y_preds <- predict(fit, as.matrix(testData[Filter(filter_out_label, names(df))]),
                       s=alpha, type="fit", mode="lambda")$fit # default mode="step"

    y_true = testData$y
    residuals = (y_true - y_preds)
    mse=sum(residuals^2)
    total_mse = total_mse + mse
  }
  if (total_mse < min_mse){
    min_mse = total_mse
    optimal_alpha = alpha
  }
}

print(paste("the optimal alpha is ", optimal_alpha))

Die Ausgabe des obigen R-Codes lautet:

> source('~.....')
[1] "the optimal alpha is  1e+06"

Unabhängig davon, was ich für die Zeile " alphas = logspace(-5, 6, 100)" eingestellt habe, erhalte ich immer den höchsten Wert für Alpha zurück.

Ich denke, hier gibt es tatsächlich zwei verschiedene Fragen:

  1. Warum ist der Alpha-Wert für die in Python geschriebene Version so instabil?

  2. Warum gibt mir die in R geschriebene Version ein anderes Ergebnis? (Mir ist klar, dass sich die logspaceFunktion von Rbis unterscheidet python, aber die eingeschriebene Version Rgibt mir immer den größten Wert alphafür den optimalen Alpha-Wert, während die Python-Version dies nicht tut).

Es wäre toll, diese Dinge zu wissen ...

Candic3
quelle
2
Ich bin mir nicht sicher, ob dies der Grund für das Problem ist, aber für das Lasso-Modell von scikit-learn (während Sie es aufrufen) müssen die Daten zentriert sein, was nicht so aussieht, als würden Sie es tun. Sie müssten den Mittelwert von x und y auf dem Trainingssatz subtrahieren und dann dieselben Werte vom Testsatz subtrahieren (zentrieren Sie die Daten nicht vor der Kreuzvalidierung oder zentrieren Sie den Testsatz mit seinem eigenen Mittelwert!). Eine Alternative besteht darin, den fit_interceptParameter beim Erstellen des Lasso-Modells zu verwenden.
user20160
Ich kann mir nicht vorstellen, dass dies die Instabilität beeinflusst, aber ich kann es versuchen ...
Candic3
3
(1) Im Python-Skript generieren Sie jedes Mal zufällige Daten, oder? Warum erwarten Sie, dass der optimale Regularisierungsparameter für alle Ihre zufälligen Ziehungen gleich ist? Unterschiedliche Daten können unterschiedliche optimale Regularisierungsparameter haben. (2) Welche Daten verwendet das R-Skript? Inwiefern unterscheiden sich die Ergebnisse des R-Skripts von Python? Sie liefern keine Ausgabe oder Vergleich.
Amöbe
1
Ich denke, die Frage ist hier kaum thematisch, da sie das Korrekturlesen von Code als wesentlichen Bestandteil der Übung beinhaltet. Wahrscheinlich sind einige der "seltsamen" Ergebnisse einfach auf Codierungsfehler zurückzuführen? Aber die Frage ist im Allgemeinen interessant. Was ist Alpha ? Zum Beispiel bin ich es gewohntλ ist die Strafintensität innerhalb von LASSO oder Ridge Regression und dann αist das Gewicht von LASSO gegenüber dem Grat in der elastischen Netzregression. Entspricht dein Alpha meinemλ?
Richard Hardy
1
Ist der Unterschied zwischen Python und R auch wirklich relevant für Ihre Hauptfrage zur Instabilität des optimalen Alphas? Indem Sie den Vergleich zwischen Python und R einbeziehen, führen Sie zusätzliche Komplexität und neue Probleme ein und maskieren damit teilweise das Wesentliche der Frage, IMHO. Der Unterschied zwischen LASSO-Implementierungen in Python und R sollte wahrscheinlich als separate Frage gestellt werden.
Richard Hardy

Antworten:

14

Ich kenne Python nicht sehr gut, aber ich habe ein Problem mit Ihrem R-Code gefunden.

Sie haben die 2 Zeilen:

residuals = sum(y_true - y_preds)
mse=residuals^2

Was die Residuen summiert, quadriert sie dann. Dies unterscheidet sich stark vom Quadrieren der Residuen und anschließenden Summieren (was anscheinend der Python-Code korrekt macht). Ich würde vermuten, dass dies ein großer Teil des Unterschieds zwischen dem R-Code und dem Python-Code ist. Korrigieren Sie den R-Code und führen Sie ihn erneut aus, um festzustellen, ob er sich eher wie der Python-Code verhält.

Ich würde auch vorschlagen, dass Sie nicht nur das "beste" Alpha und das entsprechende mse speichern, sondern alle speichern und die Beziehung zeichnen. Es kann sein, dass es für Ihr Setup eine Region gibt, die ziemlich flach ist, so dass der Unterschied zwischen der mse an verschiedenen Punkten nicht sehr groß ist. Wenn dies der Fall ist, können sehr geringfügige Änderungen an den Daten (sogar die Reihenfolge in der Kreuzvalidierung) ändern, welcher Punkt unter vielen, die im Wesentlichen gleich sind, das Minimum ergibt. Eine Situation, die zu einem flachen Bereich um das Optimum führt, führt häufig zu dem, was Sie sehen, und die Darstellung aller Alpha-Werte mit den entsprechenden mse-Werten kann aufschlussreich sein.

Greg Snow
quelle
Ihr erster Kommentar war ein toller Fang - danke. Das Problem besteht weiterhin, nachdem ich diesen Fehler behoben habe. Lassen Sie mich Ihren zweiten Vorschlag versuchen.
Candic3
@ Candic3 das ist ein toller Vorschlag. Da beide Algorithmen deterministisch sind, sollte der Startwert, wenn Sie ihn reparieren, in der Lage sein, den Lösungspfad mit dem kleinsten Winkel genau mit Ihrer DIY-Version zu reproduzieren.
Shadowtalker
Sie würden beide nur dann den gleichen Lösungspfad erzeugen, wenn die Schrittgröße genau gleich ist. Außerdem verfügt die sklearnVersion über eine integrierte Kreuzvalidierung, wie @JennyLu hervorhob, sodass ein etwas anderer Fehler auftritt.
Candic3
6

sklearn hat ein Beispiel, das fast identisch mit dem ist, was Sie hier versuchen: http://scikit-learn.org/stable/auto_examples/exercises/plot_cv_diabetes.html

In der Tat zeigt dieses Beispiel, dass Sie für jede der drei in diesem Beispiel durchgeführten Falten sehr unterschiedliche Ergebnisse für Alpha erhalten. Dies bedeutet, dass Sie der Auswahl von Alpha nicht vertrauen können, da dies eindeutig stark davon abhängt, welchen Teil Ihrer Daten Sie zum Trainieren und Auswählen von Alpha verwenden.

Ich denke nicht, dass Sie Kreuzvalidierung als etwas betrachten sollten, das "konvergiert", um Ihnen eine perfekte Antwort zu geben. Eigentlich denke ich, dass es konzeptionell fast das Gegenteil von Konvergenz ist. Sie trennen Ihre Daten und gehen für jede Falte in eine "separate Richtung". Die Tatsache, dass Sie je nach Partitionierung Ihrer Test- und Trainingsdaten unterschiedliche Ergebnisse erhalten, sollte Ihnen sagen, dass die Konvergenz zu einem perfekten Ergebnis unmöglich - und auch nicht wünschenswert ist. Der einzige Weg, um jederzeit einen konsistenten Alpha-Wert zu erhalten, besteht darin, alle Ihre Daten für das Training zu verwenden. Wenn Sie dies jedoch tun würden, würden Sie das beste Lernergebnis, aber das schlechteste Validierungsergebnis erhalten.

Jenny Lu
quelle
1
Ihr Kommentar zur x-Validierung ist interessant - ich folge nicht ganz. Ich dachte, die X-Validierung würde zur Auswahl eines Hyperparameters verwendet. Wenn die x-Validierung nicht konvergiert, was würden Sie dann zur Auswahl eines Hyperparameters verwenden?
Candic3
In dem sklearnBeispiel, das Sie auf plot_cv_diabetes zitiert haben, hat es so wenige Datenpunkte (150), dass ich nicht davon überzeugt wäre, dass dieαallein aufgrund dieses Beispiels wäre es instabil.
Candic3
Dies ist eigentlich eine nette Variante der Kreuzvalidierung. @ Candic3 Wenn es nicht konvergiert, probierst du etwas anderes aus. Das ist so, als würde ich dir sagen, dass du kein Auto über den See fahren kannst, und dann beschwerst du dich "aber ich muss rüber!" Finde eine Brücke oder gehe herum
Shadowtalker
@ Candic3 Ich habe es schnell mit allen Daten aus dem Diabetes-Datensatz (442 Punkte) ausgeführt und dies sind die Ergebnisse: [Falte 0] Alpha: 0,00010, Punktzahl: 0,50126 [Falte 1] Alpha: 0,10405, Punktzahl: 0,48495 [Falte 2] Alpha: 0,04520, Punktzahl: 0,50332
Jenny Lu
1
@JennyLu So wie ich es verstehe k folds cross-validationund (so wie ich es im obigen Beispiel gemacht habe) der Wert vonαsollte für alle Falten gleich sein. Der Wert des Parameters, den Sie schätzen (Fehler, Punktzahl oderR2, MSE usw.) sollte sich zwischen den Falten ändern. Denn versucht k folds cross-validationim Wesentlichen, den bedingten Mittelwert des von Ihnen geschätzten Parameters (den Schätzer) zu berechnen. Also denke ich nicht den Wert vonαsollte zwischen den Falten wechseln.
Candic3
5

Die Multikollinearität in x1und x2macht das ausαWert instabil im Python-Code. Die Varianz ist für die Verteilungen, die diese Variablen erzeugen, so gering, dass die Varianz der Koeffizienten aufgeblasen wird. Der Varianzinflationsfaktor (VIF) könnte berechnet werden, um dies zu veranschaulichen. Nachdem die Varianz von erhöht ist

x1 = range(n) + norm.rvs(0, 1, n) + 50
x2 =  map(lambda aval: aval*x1x2corr, x1) + norm.rvs(0, 2, n) + 500

....zu....

x1 = range(n) + norm.rvs(0, 100, n) + 50
x2 =  map(lambda aval: aval*x1x2corr, x1) + norm.rvs(0, 200, n) + 500

dann ist die α Wert stabilisiert sich.

Das Problem, dass sich der RCode vom Python-Code unterscheidet, ist jedoch immer noch ein Rätsel ...

Sother
quelle
Danke - das war's. Lassen Sie mich den Unterschied zwischen Rund Pythonetwas mehr untersuchen.
Candic3
@ Candic3, weil sie unterschiedliche Implementierungen verwenden, die bei schlecht konditionierten Problemen unterschiedliche Fehlermodi aufweisen. Wenn Sie die Dokumentation lesen, verwendet Python Lasso den Koordinatenabstieg, den Sie mit einer LAR-Lösung in R
shadowtalker vom
2

Ich werde den R-Code kommentieren:

Sie setzen Variablen an den falschen Stellen zurück, dh die Variablen min_msesollten Infaußerhalb der forSchleife optimal_alphainitialisiert und NULLdort initialisiert werden. Dies wird:

library(glmnet)
library(lars)
library(pracma)

set.seed(1)
k = 2 # number of features selected 

n = 100

x1x2corr = 1.1
x1 = seq(n) + rnorm(n, 0, 1) + 50
x2 =  x1*x1x2corr + rnorm(n, 0, 2) + 500
y = x1 + x2 +rnorm(n,0,0.5)
df = data.frame(x1 = x1, x2 = x2, y = y)
filter_out_label <- function(col) {col!="y"}

alphas = logspace(-5, 6, 50)

###
# INITIALIZE here before loop
###
min_mse = Inf
optimal_alpha = NULL
# Let's store the mse values for good measure
my_mse = c()

for (alpha in alphas){
  k = 10
  folds <- cut(seq(1, nrow(df)), breaks=k, labels=FALSE)
  # DO NOT INITIALIZE min_mse and optimal_alpha here, 
  # then you cannot find them...
  total_mse = 0
  for(i in 1:k){
    # Segement your data by fold using the which() function
    testIndexes <- which(folds==i, arr.ind=TRUE)
    testData <- df[testIndexes, ]
    trainData <- df[-testIndexes, ]

    fit <- lars(as.matrix(trainData[Filter(filter_out_label, names(df))]),
                trainData$y,
                type="lasso")
    # predict
    y_preds <- predict(fit, as.matrix(testData[Filter(filter_out_label,
                       names(df))]),
                       s=alpha, type="fit", mode="lambda")$fit 

    y_true = testData$y
    residuals = (y_true - y_preds)
    mse=sum(residuals^2)
    total_mse = total_mse + mse
  }
  # Let's store the MSE to see the effect
  my_mse <- c(my_mse, total_mse)
  if (total_mse < min_mse){
    min_mse = total_mse
    optimal_alpha = alpha
    # Let's observe the output
    print(min_mse)
  }
}

print(paste("the optimal alpha is ", optimal_alpha))
# Plot the effect of MSE with varying alphas
plot(my_mse)

Die Ausgabe sollte nun konsistent die kleinsten Alpha-Werte sein, da die Prädiktoren eine starke Kolinearität aufweisen und die Antwort nur aus den verfügbaren Prädiktoren besteht, dh es gibt keine redundante Variable, die LASSO auf Null setzen soll, in diesem Fall wir Ich möchte keine Regularisierung durchführen, dh das Kleinste alphasollte das Beste sein. Sie können die Wirkung von MSE hier sehen:

Wirkung auf mse

Beachten Sie, dass ich 50 Alphas auf derselben Skala wie Sie verwende. Bei einem Alpha-Index von 35 werden beide Variablen auf Null gesetzt, was bedeutet, dass das Modell immer dasselbe tut und die mse stagniert.

Ein besseres Problem, um MSE, CV und LASSO zu studieren

Das obige Problem ist für den LASSO nicht sehr interessant. Das LASSO führt die Modellauswahl durch, daher möchten wir sehen, dass es tatsächlich die interessierenden Parameter auswählt. Es ist beeindruckender zu sehen, dass das Modell tatsächlich ein Alpha auswählt, das die MSE tatsächlich senkt, dh uns bessere Vorhersagen gibt, indem einige Variablen weggeworfen werden. Hier ist ein besseres Beispiel, in dem ich eine Reihe redundanter Prädiktoren hinzufüge.

set.seed(1)
k = 100 # number of features selected 

n = 100

x1x2corr = 1.1
x1 = seq(n) + rnorm(n, 0, 1) + 50
x2 =  x1*x1x2corr + rnorm(n, 0, 2) + 500
# Rest of the variables are just noise
x3 = matrix(rnorm(k-2,0,(k-2)*n),n,k-2)
y = x1 + x2 +rnorm(n,0,0.5)
df = data.frame(x1 = x1, x2 = x2, y = y)
df <- cbind(df,x3)
filter_out_label <- function(col) {col!="y"}

alphas = logspace(-5, 1.5, 100)
min_mse = Inf
optimal_alpha = NULL
my_mse = c()

Dann führen Sie einfach die for-Schleife wie im obigen Code aus! Beachten Sie, dass ich das Maximum alphasvon 6 von 6 auf 1,5 gesetzt habe, um den Effekt in der folgenden Darstellung zu sehen. Jetzt ist der beste alphaWert nicht der niedrigste, aber Sie können im Diagramm sehen, dass die Kreuzvalidierungs-MSE einen Tropfen nimmt und am Ende wieder ansteigt. Der niedrigste Punkt in diesem Diagramm entspricht dem Alpha-Index mit dem niedrigsten CV-Fehler.

Besseres Lebenslaufproblem für LASSO

Gumeo
quelle