7
$\begingroup$

To explore how the LASSO regression works, I wrote a small piece of code that should optimize LASSO regression by picking the best alpha parameter.

I cannot figure out why the LASSO regression is giving me such unstable results for the alpha parameter after cross validation.

Here is my 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()

If you run this code repeatedly, it gives wildly different results for alpha:

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

Why is the alpha value not converging properly? I know that my data is synthetic, but the distribution is the same. Also, the variation is very small in x1 and x2.

what could be causing this to be so unstable?

The same thing written in R gives different results - it always returns the highest possible value for alpha as the "optimal_alpha".

I also wrote this in R, which gives me a slightly different answer, which I don't know why?

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

The output from the R code above is:

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

In fact, no matter what I set for the line "alphas = logspace(-5, 6, 100)", I always get back the highest value for alpha.

I guess there are actually 2 different questions here :

  1. Why is the alpha value so unstable for the version written in Python?

  2. Why does the version written in R give me a different result? (I realize that the logspace function is different from R to python, but the version written in R always gives me the largest value of alpha for the optimal alpha value, whereas the python version does not).

It would be great to know these things...

$\endgroup$
  • 2
    $\begingroup$ Not sure if this is what's causing the problem, but scikit-learn's lasso model (as you're invoking it) requires the data to be centered, which it doesn't look like you're doing. You'd have to subtract the mean of x and y on the training set, then subtract these same values from the test set (don't center the data before cross validation, or center the test set using its own mean!). An alternative is to use the fit_intercept parameter when you construct the lasso model. $\endgroup$ – user20160 Jul 4 '16 at 12:03
  • $\begingroup$ I can't imagine this affecting the instability, but I can try it... $\endgroup$ – Candic3 Jul 4 '16 at 20:50
  • 3
    $\begingroup$ (1) In the Python script you are generating some random data every time, right? Why do you expect that the optimal regularization parameter will be the same for all your random draws? Different data can have different optimal regularization parameters. (2) What data is the R script using? In what sense are the results from the R script different fromi Python? You don't provide any output or comparison. $\endgroup$ – amoeba says Reinstate Monica Jul 7 '16 at 11:43
  • 1
    $\begingroup$ I think the question is hardly on topic here as it involves proof-reading code as an essential part of the exercise. Probably some of the "strange" results are simply due to coding errors? But the question is interesting in general. Also, what is alpha? For example, I am used to $\lambda$ being the penalty intensity within LASSO or ridge regression and then $\alpha$ being the weight of LASSO versus ridge in elastic net regression. Does your alpha correspond to my $\lambda$? $\endgroup$ – Richard Hardy Jul 9 '16 at 12:36
  • 1
    $\begingroup$ Also, is the difference between Python and R really relevant for your main question on instability of optimal alpha? By including the comparison between Python and R you are introducing extra complexity and new problems, and thus partly masking the essence of the question, IMHO. The difference between LASSO implementations in Python and R should probably be posed as a separate question. $\endgroup$ – Richard Hardy Jul 9 '16 at 12:39
14
+50
$\begingroup$

I don't know python very well, but I did find one problem with your R code.

You have the 2 lines:

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

Which sums the residuals, then squares them. This is very different from squaring the residuals, then summing them (which it appears that the python code does correctly). I would suspect that this may be a big part of the difference between the R code and the python code. Fix the R code and run it again to see if it behaves more like the python code.

I would also suggest that instead of just saving the "best" alpha and the corresponding mse that you store all of them and plot the relationship. It could be that for your setup there is a region that is quite flat so that the difference between the mse at different points is not very big. If this is the case, then very minor changes to the data (even the order in the cross-validation) can change which point, among many that are essentially the same, gives the minimum. Having a situation that results in a flat region around the optimum will often lead to what you are seeing and the plot of all the alpha values with the corresponding mse values could be enlightening.

$\endgroup$
  • $\begingroup$ Your first comment was a great catch - thank you. The problem still persists after I fix that bug. Let me try your second suggestion. $\endgroup$ – Candic3 Jul 8 '16 at 5:06
  • $\begingroup$ @Candic3 this is a great suggestion. Also, because both algorithms are deterministic, if you fix the seed should be able to reproduce the least angle solution path exactly with your DIY version. $\endgroup$ – shadowtalker Jul 12 '16 at 12:04
  • $\begingroup$ They would both produce the same solution path only if the step size is exactly the same. Also, the sklearn version has a built-in cross validation, as @JennyLu pointed out, so it will produce an error that is a little bit different. $\endgroup$ – Candic3 Jul 18 '16 at 3:25
6
$\begingroup$

sklearn has an example that is almost identical to what you're trying to do here: http://scikit-learn.org/stable/auto_examples/exercises/plot_cv_diabetes.html

Indeed this example shows that you do get wildly varying results for alpha for each of the three folds done in that example. This means that you cannot trust the selection of alpha because it clearly is highly dependent on what portion of your data you are using to train and select alpha.

I don't think you should think of cross validation as something that will 'converge' to give you a perfect answer. Actually, I think that conceptually it is almost the opposite of converging. You are separating your data and for each fold you are going in a 'separate direction'. The fact that you get different results depending on how you partition your testing and training data should tell you that converging on one perfect result is impossible - and also not desirable. The only way you would get a consistent alpha value all the time is if you were to use all your data for training. However, if you were to do this you would get the best learning result but the worst validation result.

$\endgroup$
  • 1
    $\begingroup$ Your comment on x-validation is interesting - I don't quite follow. I thought x-validation would be used to selecting a hyperparameter. If x-validation doesn't converge, then what would you use for selecting a hyperparameter? $\endgroup$ – Candic3 Jul 9 '16 at 20:38
  • $\begingroup$ In the sklearn example you cited on plot_cv_diabetes, it has so few data points (150) that I would not be convinced that the $\alpha$'s would be unstable based on that example alone. $\endgroup$ – Candic3 Jul 9 '16 at 21:38
  • $\begingroup$ This is actually a nice take on cross validation. @Candic3 if it doesn't converge, you try something else. That's like me telling you you can't drive a car across the lake, and then you complaining "but I need to get across!" Find a bridge, or go around $\endgroup$ – shadowtalker Jul 12 '16 at 12:03
  • $\begingroup$ @Candic3 I ran it quickly with all the data from the diabetes dataset (442 points) and these are the results: [fold 0] alpha: 0.00010, score: 0.50126 [fold 1] alpha: 0.10405, score: 0.48495 [fold 2] alpha: 0.04520, score: 0.50332 $\endgroup$ – Jenny Lu Jul 13 '16 at 19:36
  • 1
    $\begingroup$ @JennyLu The way I understand k folds cross-validation, and (the way I did it in the example above) the value of $\alpha$ should be the same for all of the folds. The value of whatever parameter you're estimating (error, score or $R^2$, MSE, etc.), is what should change between the folds. Because, k folds cross-validation is essentially trying to compute the conditional mean of the parameter you're estimating (the estimand). So, I don't think the value of $\alpha$ should be changing between folds. $\endgroup$ – Candic3 Jul 18 '16 at 3:50
5
$\begingroup$

The multi-collinearity in x1 and x2 is what makes the $\alpha$ value unstable in the Python code. The variance is so small for the distributions that generate these variables, so that the variance of the coefficients are inflated. The variance inflation factor (VIF) could be computed to illustrate this. After the variance is increased from

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

....to....

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

then the $\alpha$ value stabilizes.

The issue with the R code being different from the Python code is still a mystery however...

$\endgroup$
  • $\begingroup$ Thank you - that's it. Let me explore the difference between R and Python a little more. $\endgroup$ – Candic3 Jul 10 '16 at 18:01
  • $\begingroup$ @Candic3 because they use different implementations, which have different failure modes under ill-conditioned problems. If you read the documentation, Python lasso uses coordinate descent, which you're comparing to a LAR solution in R $\endgroup$ – shadowtalker Jul 12 '16 at 12:00
2
$\begingroup$

I am going to comment on the R code:

You are resetting variables in the wrong places, i.e., the variables min_mse should be initialized as Inf outside the for loop and optimal_alpha should be initialized as NULL there. This becomes:

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)

The output should now be consistently the smallest values of alpha, because there is strong colinearity in the predictors and the response is only built from the available predictors, i.e. there are no redundant variable that we want LASSO to put to zero, in this case we do not want to perform regularization, i.e. the smallest alpha should be the best. You can see the effect of MSE here:

effect on mse

Note that I am using 50 alphas on the same scale as you. Around alpha indexed 35 both variables are slammed to zero, meaning that the model is always doing the same thing and the mse stagnates.

A better problem to study MSE, CV and the LASSO

The problem above is not very interesting for the LASSO. The LASSO performs model selection, so we want to see it actually pick out the parameters of interest. It is more impressive to see that the model is actually picking out an alpha that actually lowers the MSE, i.e. gives us better predictions by throwing out some variables. Here is a better example, where I add a bunch of redundant predictors.

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

Then you just run the for loop like in the code above! Note that I put the max of the alphas down to 1.5 from 6, just to see the effect in the plot below. Now the best alpha value is not the lowest one, but you can see in the plot that the cross-validation MSE is taking a drop and spikes up again in the end. The lowest point on that graph, corresponds to the alpha index with the lowest CV-error.

Better CV problem for LASSO

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.