Given $y \in \mathbb{R}^{n \times 1}, X \in \mathbb{R}^{n \times p}$, $p > n$, assume a LASSO-type optimization problem in the form of
$$ \hat\beta=\underset{\beta}{\operatorname{argmin}}\frac{1}{2}\left\|y - X \beta \right\|_{2}^{2} + \sum_{i}p(|\beta_i|; \gamma; \lambda)$$
where $p(|\beta_i|; \gamma; \lambda)$ is a non-convex function in $\beta$ with $\gamma$ and $\lambda$ denoting the degree of regularization and non-convexity. Typically known examples are SCAD, MCP penalties, among probably many others.
My question: I want to implement solutions for known penalties and test out some custom penalties. What are currently the faster known solvers that would be able to adequately tackle such problems? Which family of solvers should I be looking at for quick/feasible computations? Are there any general solvers, or do they highly depend on the specific penalty function?
My initial approach was to start from Augmented Lagrangian for every $\lambda, \gamma$ throughout the solution path, which seem to work for small $n$ and $p$ values, given $p < n$, but doesn't seem to scale well into bigger problems.
