def gradient(self):
        r"""Gradient operator of the functional.

        The gradient of the Huber functional is given by

            .. math::
                \nabla f_{\gamma}(x) =
                \begin{cases}
                \frac{1}{\gamma} x & \text{if } \|x\|_2 \leq \gamma \\
                \frac{1}{\|x\|_2} x & \text{else}
                \end{cases}.

        Examples
        --------
        Check that the gradient norm is less than the norm of the one element:

        >>> space = odl.uniform_discr(0, 1, 14)
        >>> norm_one = space.one().norm()
        >>> x = odl.phantom.white_noise(space)
        >>> huber_norm = odl.solvers.Huber(space, gamma=0.1)
        >>> grad = huber_norm.gradient(x)
        >>> tol = 1e-5
        >>> grad.norm() <=  norm_one + tol
        True

        Redo previous example for a product space in two dimensions:

        >>> domain = odl.uniform_discr([0, 0], [1, 1], [5, 5])
        >>> space = odl.ProductSpace(domain, 2)
        >>> norm_one = space.one().norm()
        >>> x = odl.phantom.white_noise(space)
        >>> huber_norm = odl.solvers.Huber(space, gamma=0.2)
        >>> grad = huber_norm.gradient(x)
        >>> tol = 1e-5
        >>> grad.norm() <=  norm_one + tol
        True
        """

        functional = self

        class HuberGradient(Operator):

            """The gradient operator of this functional."""

            def __init__(self):
                """Initialize a new instance."""
                super(HuberGradient, self).__init__(
                    functional.domain, functional.domain, linear=False)

            def _call(self, x):
                """Apply the gradient operator to the given point."""
                if isinstance(self.domain, ProductSpace):
                    norm = PointwiseNorm(self.domain, 2)(x)
                else:
                    norm = x.ufuncs.absolute()

                grad = x / functional.gamma

                index = norm.ufuncs.greater_equal(functional.gamma)
                if isinstance(self.domain, ProductSpace):
                    for xi, gi in zip(x, grad):
                        gi[index] = xi[index] / norm[index]
                else:
                    grad[index] = x[index] / norm[index]

                return grad

        return HuberGradient()