  This paper concerns with the group zero-norm regularized least squares
estimator which, in terms of the variational characterization of the zero-norm,
can be obtained from a mathematical program with equilibrium constraints
(MPEC). By developing the global exact penalty for the MPEC, this estimator is
shown to arise from an exact penalization problem that not only has a favorable
bilinear structure but also implies a recipe to deliver equivalent DC
estimators such as the SCAD and MCP estimators. We propose a multi-stage convex
relaxation approach (GEP-MSCRA) for computing this estimator, and under a
restricted strong convexity assumption on the design matrix, establish its
theoretical guarantees which include the decreasing of the error bounds for the
iterates to the true coefficient vector and the coincidence of the iterates
after finite steps with the oracle estimator. Finally, we implement the
GEP-MSCRA with the subproblems solved by a semismooth Newton augmented
Lagrangian method (ALM) and compare its performance with that of SLEP and
MALSAR, the solvers for the weighted $\ell_{2,1}$-norm regularized estimator,
on synthetic group sparse regression problems and real multi-task learning
problems. Numerical comparison indicates that the GEP-MSCRA has significant
advantage in reducing error and achieving better sparsity than the SLEP and the
MALSAR do.
