1
$\begingroup$

Summary

When I run glmnet on the same LASSO problem by successively decreasing the convergence tolerance (threshold), I have observed that in some cases, the objective function values increase as the tolerance decreases. This does not seem right given that coordinate descent, which is the algorithm glmnet uses to solve the LASSO, is guaranteed to decrease the objective function value at every iteration. This is mentioned, for example, in Wright's manuscript.

Details

The Lasso problem is defined as follows:

$$ \min_{\beta \in \mathbb{R}^p} \frac{1}{2n} \Vert y- X\beta \Vert _2^2 + \lambda\Vert\beta \Vert_1, \; \mathrm{where} \; X \in \mathbb{R}^{n \times p}, y \in \mathbb{R}^n \; \mathrm{and} \; \lambda > 0 $$

glmnet gives you the solution to the above problem for a sequence of $\lambda_i$ values. The R code below will first generate a small problem and then successively run glmnet on the same problem with decreased tolerance values. It then takes the $\tilde{\beta^i}$s, the solutions corresponding to $\lambda_i$, and computes the objective function value. It performs this for every glmnet run and then computes the difference in the objective. Surprisingly, there are quite a few $\lambda_i$s that yield a higher objective function value despite the lower tolerance value used during the computation.

# Load libraries
library(Matrix)
library(glmnet)

# Generate a small problem
# Problem parameters
num_rows = 5
num_cols = 19
sparsity = 0.2
# Set the seed so runs stay the same
set.seed(num_cols)
# Create the X matrix
X = Matrix(rnorm(num_rows*num_cols), nrow=num_rows, ncol=num_cols)
# Create an explanatory vector b that is used to generate y
b = Matrix(0,nrow=ncol(X),ncol=1)
# Set the seed so runs stay the same
set.seed(num_cols)
# Determine nonzero vars
nzvars = sample.int(num_cols, size = max(num_rows, ceiling(sparsity*num_cols)) , replace = FALSE, prob = NULL)
b[nzvars,1] = 1
b = Matrix(b)
# Create the true response vector y
y = X %*% b
# Set the seed so runs stay the same
set.seed(num_cols)
# Add noise to the response vector y
y = apply(y, 2, jitter)
y = Matrix(y)


# Test glmnet 
tol=1.0e-3
while(tol >= 1.0e-5) {
  
  # Run glmnet with tol
  cat(paste("Running glmnet with tolerance: ", tol, "\n"))
  res <- glmnet(X, y, thresh = tol)
  # Compute objective function values for each lambda
  obj_val1 <- numeric(ncol(res$beta))
  for(i in 1:ncol(res$beta)) {
    obj_val1[i] = (1.0/(2.0*nrow(X)))*crossprod(y - X %*% as.matrix(res$beta[,i])) + res$lambda[i] * sum(abs(res$beta[,i]))
  }
  
  # Run glmnet again this time with tol/10.0
  tol = tol/10.0
  cat(paste("Running glmnet with tolerance: ", tol, "\n"))
  res <- glmnet(X, y, thresh = tol)
  # Compute objective function values for each lambda
  obj_val2 <- numeric(ncol(res$beta))
  for(i in 1:ncol(res$beta)) {
    obj_val2[i] = (1.0/(2.0*nrow(X)))*crossprod(y - X %*% as.matrix(res$beta[,i])) + res$lambda[i] * sum(abs(res$beta[,i]))
  }
  
  # Compute the difference between the two objectives 
  dim = min(length(obj_val1), length(obj_val2))
  obj_val_diff = obj_val1[1:dim] - obj_val2[1:dim]
  df <- data.frame(
    lambda = res$lambda[which(obj_val_diff < 0)],
    obj_val_diff = obj_val_diff[which(obj_val_diff < 0)],
    tol1 = rep(10.0*tol, length(which(obj_val_diff < 0))),
    obj_val1 = obj_val1[which(obj_val_diff < 0)],
    tol2 = rep(tol, length(which(obj_val_diff < 0))),
    obj_val2 = obj_val2[which(obj_val_diff < 0)]
  )
  print(df)
}

For tolerances $10^{-3}$ and $10^{-4}$ there is no increase in the objective. However, when going from $10^{-4}$ to $10^{-5}$ or from $10^{-5}$ to $10^{-6}$ there is a clear increase in the objective for some of the lambdas. Below is the output for $10^{-5}$ to $10^{-6}$:

lambda      obj_val_diff    tol1    obj_val1    tol2    obj_val2
0.7892003   -2.157485e-04   1e-05   1.636470    1e-06   1.636686
0.7533300   -3.773846e-05   1e-05   1.645333    1e-06   1.645371
0.7190900   -1.523336e-05   1e-05   1.651397    1e-06   1.651412
0.5970011   -6.661338e-16   1e-05   1.654829    1e-06   1.654829
0.5698665   -1.686165e-04   1e-05   1.653610    1e-06   1.653778
0.5439652   -1.610012e-04   1e-05   1.650942    1e-06   1.651103
0.5192411   -1.565519e-04   1e-05   1.647132    1e-06   1.647288
0.4956408   -1.520903e-04   1e-05   1.642348    1e-06   1.642500
0.4731132   -1.475959e-04   1e-05   1.636737    1e-06   1.636885
0.4516094   -1.430910e-04   1e-05   1.630430    1e-06   1.630573
0.4310831   -1.385951e-04   1e-05   1.623541    1e-06   1.623680
0.4114897   -1.312556e-04   1e-05   1.616189    1e-06   1.616321
0.3927869   -1.463182e-04   1e-05   1.607602    1e-06   1.607748

Question

I have no explanation for the above occurrence. Does anybody know why this happens?

$\endgroup$

1 Answer 1

0
$\begingroup$

I am working on it right now. My ADMM results for lasso does not converge to glmnet results. After a couple of frustration days, I have noticed that glmnet might be using different objective function compared to the one you gave. Seems like they might be dividing L2 norm by number of observations (see function description).

$\endgroup$
1
  • $\begingroup$ Hello, both in the formulation and the code above I'm dividing by the number of observations. (number of observations = n = nrow(X)). GLMNET's termination criteria is an objective function value based termination so it's very flaky, whereas ADMM uses the sub-gradient, which gives a more accurate result. To get GLMNET aligned with traditional gradient/sub-gradient based methods you need to crank up the thresh (i.e. tolerance) to something absurd like 1e-15. $\endgroup$
    – gomfy
    Sep 28, 2020 at 11:14

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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