def derivative(self, x):
        """Derivative of the product space operator.

        Parameters
        ----------
        x : `domain` element
            The point to take the derivative in

        Returns
        -------
        adjoint : linear`ProductSpaceOperator`
            The derivative

        Examples
        --------
        >>> r3 = odl.rn(3)
        >>> pspace = odl.ProductSpace(r3, r3)
        >>> I = odl.IdentityOperator(r3)
        >>> x = pspace.element([[1, 2, 3], [4, 5, 6]])

        Example with linear operator (derivative is itself)

        >>> prod_op = ProductSpaceOperator([[0, I], [0, 0]],
        ...                                domain=pspace, range=pspace)
        >>> prod_op(x)
        ProductSpace(rn(3), 2).element([
            [ 4.,  5.,  6.],
            [ 0.,  0.,  0.]
        ])
        >>> prod_op.derivative(x)(x)
        ProductSpace(rn(3), 2).element([
            [ 4.,  5.,  6.],
            [ 0.,  0.,  0.]
        ])

        Example with affine operator

        >>> residual_op = I - r3.element([1, 1, 1])
        >>> op = ProductSpaceOperator([[0, residual_op], [0, 0]],
        ...                           domain=pspace, range=pspace)

        Calling operator gives offset by [1, 1, 1]

        >>> op(x)
        ProductSpace(rn(3), 2).element([
            [ 3.,  4.,  5.],
            [ 0.,  0.,  0.]
        ])

        Derivative of affine operator does not have this offset

        >>> op.derivative(x)(x)
        ProductSpace(rn(3), 2).element([
            [ 4.,  5.,  6.],
            [ 0.,  0.,  0.]
        ])
        """
        # Lazy import to improve `import odl` time
        import scipy.sparse

        # Short circuit optimization
        if self.is_linear:
            return self

        deriv_ops = [op.derivative(x[col]) for op, col in zip(self.ops.data,
                                                              self.ops.col)]
        data = np.empty(len(deriv_ops), dtype=object)
        data[:] = deriv_ops
        indices = [self.ops.row, self.ops.col]
        shape = self.ops.shape
        deriv_matrix = scipy.sparse.coo_matrix((data, indices), shape)
        return ProductSpaceOperator(deriv_matrix, self.domain, self.range)