def proximal_arg_scaling(prox_factory, scaling):
    r"""Calculate the proximal of function F(x * scaling).

    Parameters
    ----------
    prox_factory : callable
        A factory function that, when called with a step size, returns the
        proximal operator of ``F``
    scaling : float or sequence of floats or space element
        Scaling parameter. The permissible types depent on the stepsizes
        accepted by prox_factory. It may not contain any nonzero imaginary
        parts. If it is a scalar, it may be zero, in which case the
        resulting proxmial operator is the identity. If not a scalar,
        it may not contain any zero components.

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

    Notes
    -----
    Given a functional :math:`F`, and scaling factor :math:`\alpha` this is
    calculated according to the rule

    .. math::
        \mathrm{prox}_{\sigma F(\alpha \, \cdot)}(x) =
        \frac{1}{\alpha}
        \mathrm{prox}_{\sigma \alpha^2 F(\cdot) }(\alpha x)

    where :math:`\sigma` is the step size.

    For reference on the identity used, see [CP2011c].

    References
    ----------
    [CP2011c] Combettes, P L, and Pesquet, J-C. *Proximal splitting
    methods in signal processing.* In:  Bauschke, H H, Burachik, R S,
    Combettes, P L, Elser, V, Luke, D R, and Wolkowicz, H. Fixed-point
    algorithms for inverse problems in science and engineering, Springer,
    2011.
    """

    # To begin, we could check for two things:
    # * Currently, we do not support complex scaling. We could therefore catch
    #   nonempty imaginary parts.
    # * If some components of scaling are zero, then the following routine will
    #   crash with a division-by-zero error. The correct solution would be to
    #   just keep these components and do the following computations only for
    #   the others.
    # Since these checks are computationally expensive, we do not execute them
    # unconditionally, but only if the scaling factor is a scalar:
    if np.isscalar(scaling):
        if scaling == 0:
            return proximal_const_func(prox_factory(1.0).domain)
        elif scaling.imag != 0:
            raise ValueError("Complex scaling not supported.")
        else:
            scaling = float(scaling.real)
    else:
        scaling = np.asarray(scaling)

    def arg_scaling_prox_factory(sigma):
        """Create proximal for the translation with a given sigma.

        Parameters
        ----------
        sigma : positive float
            Step size parameter

        Returns
        -------
        proximal : `Operator`
            The proximal operator of ``sigma * F( . * a)`` where ``sigma`` is
            the step size
        """
        scaling_square = scaling * scaling
        prox = prox_factory(sigma * scaling_square)
        space = prox.domain
        mult_inner = MultiplyOperator(scaling, domain=space, range=space)
        mult_outer = MultiplyOperator(1 / scaling, domain=space, range=space)
        return mult_outer * prox * mult_inner

    return arg_scaling_prox_factory