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

• 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. Sep 28, 2020 at 11:14