# Nonzero coefficient ordering in least-square LASSO

Consider least-square LASSO over standardized training data $$(\boldsymbol{X},\boldsymbol{y})$$. Assume $$|\boldsymbol{x}_j\cdot\boldsymbol{y}|>|\boldsymbol{x}_k\cdot\boldsymbol{y}|$$. In other words, $$\boldsymbol{x}_j$$ is a more positively/negatively correlated column with $$\boldsymbol{y}$$ than $$\boldsymbol{x}_k$$.

Would coefficients $$\beta_j=0\implies\beta_k=0$$ ?

If yes, does it apply to least-square Elastic net?

• Because this condition on the correlations tells us almost nothing about the multivariate relationships among $(X,y),$ the answers must be no and no. Look for simple counterexamples in the form $y = x_1+x_2$ and $x_3 = y +\varepsilon$ for low-variance random values $\varepsilon.$ Using $x_1$ and $x_2$ gives a perfect fit, whence $x_3$ is superfluous; but each of those variables may be considerably less correlated with $y$ than $x_3$ is.
– whuber
Commented Feb 10, 2023 at 21:54

As @whuber commented, the answer is no and no. R demo:

set.seed(789)
Nrow = 1000
Ncol = 10
X = apply(matrix(rnorm(Nrow * Ncol), nrow = Nrow), 2, function(x)
{
x = x - mean(x)
x / sqrt(mean(x ^ 2))
})
rg = sort(c(0.2, 0.3) * sample(c(-1, 1), 1))
beta = runif(5, rg[1], rg[2])
y = X[,1:5] %*% beta # Make the first 5 columns the true predictors

# Make the last 5 columns increasingly more correlated with y.
r = 0.6
X[, 6:10] = apply(X[, 6:10], 2, function(x)
{
x = sort(x, decreasing = sample(c(T, F), 1))
ind = 1:as.integer(round(nrow(X) * r))
x[ind] = sample(x[ind])
r <<- r * 0.85; x
})
ind = order(y)
X[, 1:5] = X[ind, 1:5]
y = y[ind]
ind = order(apply(X, 2, function(x) abs(sum(x * y))))
beta = beta[ind[1:5]]
X = X[, ind]

# At this point, the first 5 columns of X are the true predictors,
# but each of them has low abs correlation with y. The last 5 columns
# are not true predictors but have high abs correlations with y. Demo:
cat("Correlation between X's columns and y: ", round(as.numeric(cor(X, y)), 3), "\n")
cat("Coefficient of determination of the first 5 columns =",
1 - mean((X[, 1:5] %*% beta - y) ^ 2) / mean(y ^ 2))

mdl = glmnet::glmnet(X, y, nlambda = 20, lambda.min.ratio = 0.01)
round(mdl\$beta, 5)