Parameter choice rules for L1 regularization?

Question

I am solving an L1 regularized least squares (LASSO) of the form:

$$ \arg \min_{\boldsymbol{x}} \frac{1}{2} {\left\| A \boldsymbol{x} - \boldsymbol{y} \right\|}_{2}^{2} + \lambda {\left\| \boldsymbol{x} \right\|}_{1} $$

I saw that for the L2 norm case there are several methods to obtain the magnitude of $\lambda$ as seen in the study Comparing parameter choice methods for regularization of ill-posed problems. For the L1 case I was told to look at AIC\BIC but I have a feeling I am not using them correctly.

For example, y is a 1x1000 vector and A is a 1000x1000 dictionary, and I plan it so that x will be a sparse vector of a length 1000 with only 5-15 non zeros values (no noise). By trial and error the range of $\lambda$s that retrieve accurate solutions is broad in the range 1e-4-1e-2 (using CVX). I want to find the optimal value for $\lambda$, similar to generalized cross validation that is done in L2 norm case (for example).

I am calculating AIC\BIC (using aicbic function in matlab) as follows: I scan (for-loop) values of $\lambda$ in the range 1e-7 to 1. For each $\lambda$ value in the loop I use the log of the optimal value of the objective function found by CVX (cvx_optval) as the log-likelihood, for the number of estimated model parameters I use the number of non-zero (or |x|>1e-6) entries in the solution of x, and for the sample size I simply use the length of y.

I get the minimum of the AIC\BIC estimate for a $\lambda$ value around 1, this is similar to the estimated maximal value of $\lambda$ given by $\lambda_{max} = norm(2*(At*y),inf)$ that is seen in Boyd's l1_ls guide. This value is far from optimal and doesn't give the correct solution (not accurate in the "position" or the # of non-zero entries expected). The reason is because the number of non-zero entries that is used for the estimated model parameters is dropping to 1-3, so the total function gets a minima there.

What am I doing wrong, and should I use AIC\BIC or something else and how?

Answer 1

For example, y is a 1x1000 vector and A is a 1000x1000 dictionary, and I plan it so that x will be a sparse vector of a length 1000 with only 5-15 non zeros values (no noise).

When tunning the $\lambda$ hyperparameter, you are gaining at minimizing the overall loss. This does not necessarily have anything to do with the number of non-zero values.

I am calculating AIC\BIC (using aicbic function in matlab) as follows: I scan (for-loop) values of $\lambda$ in the range 1e-7 to 1. For each $\lambda$ value in the loop I use the log of the optimal value of the objective function [...]

Grid search is far from a great optimization algorithm. There are infinitely many possible values in the range 1e-7 to 1, you are arbitrarily picking a small, finite number of them to loop over. You should use instead one of many optimization algorithms available (e.g. Bayesian optimization).

I get the minimum of the AIC\BIC estimate for a $\lambda$ value around 1, this is similar to the estimated maximal value of $\lambda$ given by $\lambda_{max} = norm(2*(At*y),inf)$ that is seen in [Boyd's l1_ls][3] guide. This value is far from optimal and doesn't give the correct solution (not accurate in the "position" or the # of non-zero entries expected). The reason is because the number of non-zero entries that is used for the estimated model parameters is dropping to 1-3, so the total function gets a minima there.

As stated earlier, it does what you asked it for. You asked it for giving you the value that has minima at the objective function, not to give you the specific number of the non-zero values. If you want exactly 5-15 non-zero values, your objective function should reflect that (e.g. the loss should be infinite when the number is not in the 5-15 range), but it would be a different optimization problem, and likely an ill-posed one.