  Hierarchical models with gamma hyperpriors provide a flexible,
sparse-promoting framework to bridge $L^1$ and $L^2$ regularizations in
Bayesian formulations to inverse problems. Despite the Bayesian motivation for
these models, existing methodologies are limited to \textit{maximum a
posteriori} estimation. The potential to perform uncertainty quantification has
not yet been realized. This paper introduces a variational iterative
alternating scheme for hierarchical inverse problems with gamma hyperpriors.
The proposed variational inference approach yields accurate reconstruction,
provides meaningful uncertainty quantification, and is easy to implement. In
addition, it lends itself naturally to conduct model selection for the choice
of hyperparameters. We illustrate the performance of our methodology in several
computed examples, including a deconvolution problem and sparse identification
of dynamical systems from time series data.
