  In the context of linear inverse problems, we propose and study a general
iterative regularization method allowing to consider large classes of
regularizers and data-fit terms. The algorithm we propose is based on a
primal-dual diagonal {descent} method. Our analysis establishes convergence as
well as stability results. Theoretical findings are complemented with numerical
experiments showing state of the art performances.
