def proximal_l2_squared(space, lam=1, g=None):
    r"""Proximal operator factory of the squared l2-norm/distance.

    Function for the proximal operator of the convex conjugate of the
    functional F where F is the l2-norm (or distance to g, if given)::

        F(x) =  lam ||x - g||_2^2

    with x and g elements in ``space``, scaling factor lam, and given data g.

    Parameters
    ----------
    space : `LinearSpace`
        Domain of F(x). Needs to be a Hilbert space.
        That is, have an inner product (`LinearSpace.inner`).
    lam : positive float, optional
        Scaling factor or regularization parameter.
    g : ``space`` element, optional
        An element in ``space``. Default: ``space.zero``.

    Returns
    -------
    prox_factory : function
        Factory for the proximal operator to be initialized

    Notes
    -----
    The squared :math:`L_2`-norm/distance :math:`F` is given by

    .. math::
        F(x) =  \lambda \|x - g\|_2^2.

    For a step size :math:`\sigma`, the proximal operator of :math:`\sigma F`
    is given by

    .. math::
        \mathrm{prox}_{\sigma F}(x) = \frac{x + 2 \sigma \lambda g}
        {1 + 2 \sigma \lambda}.

    See Also
    --------
    proximal_l2 : proximal without square
    proximal_convex_conj_l2_squared : proximal for convex conjugate
    """
    class ProximalL2Squared(Operator):

        """Proximal operator of the squared l2-norm/dist."""

        def __init__(self, sigma):
            """Initialize a new instance.

            Parameters
            ----------
            sigma : positive float or pointwise positive space.element
                Step size parameter. If scalar, it contains a global stepsize,
                otherwise the space.element defines a stepsize for each point.
            """
            super(ProximalL2Squared, self).__init__(
                domain=space, range=space, linear=g is None)
            if np.isscalar(sigma):
                self.sigma = float(sigma)
            else:
                self.sigma = space.element(sigma)

        def _call(self, x, out):
            """Apply the operator to ``x`` and store the result in ``out``"""
            # (x + 2*sig*lam*g) / (1 + 2*sig*lam))
            sig = self.sigma
            if np.isscalar(sig):
                if g is None:
                    out.lincomb(1 / (1 + 2 * sig * lam), x)
                else:
                    out.lincomb(1 / (1 + 2 * sig * lam), x,
                                2 * sig * lam / (1 + 2 * sig * lam), g)
            else:   # sig in space
                if g is None:
                    x.divide(1 + 2 * sig * lam, out=out)
                else:
                    if x is out:
                        # Can't write to `out` since old `x` is still needed
                        tmp = sig.multiply(2 * lam * g)
                        out.lincomb(1, x, 1, tmp)
                    else:
                        sig.multiply(2 * lam * g, out=out)
                        out.lincomb(1, x, 1, out)
                    out.divide(1 + 2 * sig * lam, out=out)

    return ProximalL2Squared