I am reading a paper which says that LASSO is asymptotically biased while SCAD is not. I take asymptotic (un)biasedness to concern the slope estimators from LASSO and SCAD as the sample size goes to infinity, but I am not sure. I wonder what exactly these statements mean, under which assumptions they hold and whether these assumptions are realistic.
In my understanding, LASSO is asymptotically biased given fixed regularization intensity $\lambda$. However, realistically $\lambda$ would not be kept fixed as the sample size grows but would rather be reduced (e.g. this would be the case if one used LOOCV for selecting $\lambda$, something that is fairly common), reducing the bias accordingly. Taking this to the limit, it appears LASSO would not be asymptotically biased.
On the other hand, if we look at some typical pictures illustrating LASSO and SCAD estimators, they often consider them as functions of the slope coefficient. There, we see that LASSO is asympt. biased while SCAD is asympt. unbiased when the slope coefficient (rather than the sample size) goes to infinity. (See below.) Again, I think $\lambda$ is fixed here.
So I am confused, and hence my question.
Update 1: Here is a small simulation compatible with the statement that LASSO is asymptotically biased. $\lambda$ is selected via 8-fold CV.
(A simulation where $\lambda$ is selected by LOOCV yields similar results for sample sizes up to $2^{15}$; larger sample sizes are not feasible on my laptop.)
The figure corresponds to Setup no. 2 in the code: one relevant regressor, one irrelevant regressor. The top figure is the bias of the coefficient on the relevant regressor, the bottom figure is the bias of the coefficient on the irrelevant regressor. As you can see, the bias of the relevant regressor shrinks initially with the sample size $n$ but stops and remains nonzero after $n$ gets sufficiently large. (It may be a bit difficult to see from the graph, but the bias indeed stops shrinking for the first regressor after $n=2^{13}$.)
library(glmnet)
m=1e3 # number of simulation runs
ns=2^c(5:17) # sample sizes
k=length(ns)
bias1=rep(NA,k); names(bias1)=ns; bias2=bias1
for(j in 1:k){
n=ns[j]; print(paste0(Sys.time()," n = ",n))
beta1_hat=beta2_hat=rep(NA,m)
for(i in 1:m){
set.seed(i); data=matrix(rnorm(3*n),ncol=3); x1=data[,1]; x2=data[,2]; eps=data[,3]
# Choose one of the 3 lines below:
#beta1=0; beta2=0 # Setup no. 1: both regressors irrelevant
beta1=1; beta2=0 # Setup no. 2: only the first regressor relevant
#beta1=1; beta2=1 # Setup no. 3: both regressors relevant
y=beta1*x1+beta2*x2+eps # zero mean for simplicity
cvfit=cv.glmnet(x=cbind(x1,x2),y=y,nfolds=8)
coef=coef(cvfit,s="lambda.min")
beta1_hat[i]=coef[2]
beta2_hat[i]=coef[3]
}
bias1[j]=mean(beta1_hat)-beta1
bias2[j]=mean(beta2_hat)-beta2
}; print(paste0(Sys.time()," The end"))
par(mfrow=c(2,1),mar=c(4,4,2,0.5))
at=c(1:k); labels=paste0("2^",log(ns,2))
ylim=range(bias1,-bias1); plot(bias1,ylim=ylim,xaxt="n",xlab="sample size"); axis(side=1,at=at,labels=labels); abline(h=0)
ylim=range(bias2,-bias2); plot(bias2,ylim=ylim,xaxt="n",xlab="sample size"); axis(side=1,at=at,labels=labels); abline(h=0)
par(mfrow=c(1,1))
Update 2: searching for "asymptotically biased" online yields Javanmard & Montanari "Confidence intervals and hypothesis testing for high-dimensional regression" (2014) where the matter is discussed in the first 12 pages of the document, especially Theorems 6-8 and some discussion that follows in Section 2.2. The material is unfortunately highly technical and presented using fairly complicated notation. Figuring it out by myself is quite a challenge.