I want to deconvolve an image $h$ by a kernel $f$. More precisely, let
$$G = \operatorname*{argmin}_g \|f \ast g - h\|_2$$
be the set of least-squares solutions. I want to find the least-norm solution in this set:
$$ \hat{g} \in \operatorname*{argmin}_{g \in G} \|g\|_2 $$
$$ \hat{g} = (f \mathrel{\ast})^+ h $$
where $^+$ is the pseudoinverse and $f \mathrel{\ast}$ is the Toeplitz matrix constructed from $f$.
What algorithm can efficiently compute $\hat{g}$, particularly in the setting of 2D images? Richardson–Lucy deconvolution doesn't work because I get division-by-zero errors (e.g. when $f$ is the discrete Laplace operator). Are there other, possibly gradient-based, iterative methods that don't suffer from this problem?