  Choosing an appropriate regularization term is necessary to obtain a
meaningful solution to an ill-posed linear inverse problem contaminated with
measurement errors or noise. The $\ell_p$ norm covers a wide range of choices
for the regularization term since its behavior critically depends on the choice
of $p$ and since it can easily be combined with a suitable regularization
matrix. We develop an efficient algorithm that simultaneously determines the
regularization parameter and corresponding $\ell_p$ regularized solution such
that the discrepancy principle is satisfied. We project the problem on a
low-dimensional Generalized Krylov subspace and compute the Newton direction
for this much smaller problem. We illustrate some interesting properties of the
algorithm and compare its performance with other state-of-the-art approaches
using a number of numerical experiments, with a special focus of the sparsity
inducing $\ell_1$ norm and edge-preserving total variation regularization.
