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?
