gradient descent method. As a matter of fact, results with non-convex regularization are difficult to reproduce because different numerical optimization procedures can lead to different local minima. Therefore the quality of the solution heavily depend on the numerical procedure used.

In the following, we consider a specific numerical procedure referred to as multi-stage convex relaxation in [@Zhang09-multistage]. The algorithm is given in Figure \[fig:multi-stage-sparse\]. The procedure converges to a local optimal solution of (\[eq:sparse-reg\]) due to a simple concave duality argument, where (\[eq:sparse-reg\]) is rewritten as $$\hat{{{\mathbf w}}} = \arg\min_{{{\mathbf w}}} \min_{\{\lambda_j \geq 0\}}\left[ \frac{1}{n}\| X {{\mathbf w}}- {{\mathbf y}}\|_2^2 + 
\sum_{j=1}^p \lambda_j |{{\mathbf w}}_j| + \sum_{j=1}^p g^*(\lambda_j) \right] ,$$ with $g^*(\lambda_j)=\max((\lambda-\lambda_j)\theta,0)$. The procedure of Figure \[fig:multi-stage-sparse\] can be regarded as an alternating optimization method to solve this joint optimization problem of ${{\mathbf w}}$ and $\{\lambda_j\}$, where the first step solves for ${{\mathbf w}}$ with $\{\lambda_j\}$ fixed, and the second step is the closed form solution of $\{\lambda_j\}$ with ${{\mathbf w}}$ fixed. A more detailed discussion can be found in [@Zhang09-multistage]. Our goal is to show that this procedure can achieve unbiased feature selection as described in [@Zhang10-mc+].

Initialize $\lambda_j^{(0)}=\lambda$ for $j=1,\ldots,d$\
For $\ell=1,2,\ldots$

-   Let $$\hat{{{\mathbf w}}}^{(\ell)} = \arg\min_{{{\mathbf 