def prox_l0(v, alpha):
    r"""Compute the proximal operator of the :math:`\ell_0` "norm" (hard
    thresholding)

     .. math::
      \mathrm{prox}_{\alpha f}(v) = \mathcal{S}_{0,\alpha}(\mathbf{v})
      = \left\{ \begin{array}{ccc} 0 & \text{if} &
      | v | < \sqrt{2 \alpha} \\ v &\text{if}  &
      | v | \geq \sqrt{2 \alpha} \end{array} \right. \;,

    where :math:`f(\mathbf{x}) = \|\mathbf{x}\|_0`. The approach taken
    here is to start with the definition of the :math:`\ell_0` "norm"
    and derive the corresponding proximal operator. Note, however, that
    some authors (e.g. see Sec. 2.3 of :cite:`kowalski-2014-thresholding`)
    start by defining the hard thresholding rule and then derive the
    corresponding penalty function, which leads to a simpler form for
    the thresholding rule and a more complicated form for the penalty
    function.

    Unlike the corresponding :func:`norm_l0`, there is no need for an
    `axis` parameter since the proximal operator of the :math:`\ell_0`
    norm is the same when taken independently over each element, or
    over their sum.


    Parameters
    ----------
    v : array_like
      Input array :math:`\mathbf{v}`
    alpha : float or array_like
      Parameter :math:`\alpha`

    Returns
    -------
    x : ndarray
      Output array
    """

    return (np.abs(v) >= np.sqrt(2.0 * alpha)) * v
r"""Compute the proximal operator of the :math:`\ell_0` "norm" (hard
    thresholding)

     .. math::
      \mathrm{prox}_{\alpha f}(v) = \mathcal{S}_{0,\alpha}(\mathbf{v})
      = \left\{ \begin{array}{ccc} 0 & \text{if} &
      | v | < \sqrt{2 \alpha} \\ v &\text{if}  &
      | v | \geq \sqrt{2 \alpha} \end{array} \right. \;,

    where :math:`f(\mathbf{x}) = \|\mathbf{x}\|_0`. The approach taken
    here is to start with the definition of the :math:`\ell_0` "norm"
    and derive the corresponding proximal operator. Note, however, that
    some authors (e.g. see Sec. 2.3 of :cite:`kowalski-2014-thresholding`)
    start by defining the hard thresholding rule and then derive the
    corresponding penalty function, which leads to a simpler form for
    the thresholding rule and a more complicated form for the penalty
    function.

    Unlike the corresponding :func:`norm_l0`, there is no need for an
    `axis` parameter since the proximal operator of the :math:`\ell_0`
    norm is the same when taken independently over each element, or
    over their sum.


    Parameters
    ----------
    v : array_like
      Input array :math:`\mathbf{v}`
    alpha : float or array_like
      Parameter :math:`\alpha`

    Returns
    -------
    x : ndarray
      Output array