Sparse signal recovery using a signal dictionary is commonly framed as a regularized least-squares problem, where sparsity is promoted through an \(\ell _1\) -norm on the coefficients. However, non-convex regularization, such as the smoothly clipped absolute deviation (SCAD) regularization, have been shown to achieve superior sparse recovery performance. To enhance robustness against outliers and heavy-tailed noise, we develop a sparse recovery model that integrates a Huber loss with the SCAD regularization. Due to the non-convexity and nonsmoothness of the problem, we propose DC programming by decomposing the SCAD regularization into the difference of convex (DC) functions. We then employ a proximal majorization-minimization (PMM) framework to handle the DC term, followed by a semismooth Newton (SSN) method to solve the resulting convex subproblem. We show that the SSN method achieves a fast local convergence rate to the subproblem under certain assumptions. To evaluate the numerical performance of PMM-SSN, we also implement two variants of the alternating direction methods of multipliers (ADMM) for comparison. Finally, we perform a series of numerical experiments using both simulated and real data to demonstrate the benefits of the regularized model and the efficiency of PMM-SSN. The results indicate that PMM-SSN is much faster than both DCA-ADMM and DCA-DADMM, while maintaining a comparable recovery accuracy.