Added L0 proximal operator

Author: mrava87

pyproximal/projection/L1.py

@@ -40,4 +40,8 @@ def __init__(self, n, radius, maxiter=100, xtol=1e-5):
         self.simplex = SimplexProj(n, radius, maxiter, xtol)
 
     def __call__(self, x):
-        return np.sign(x) * self.simplex(np.abs(x))
+        if np.iscomplexobj(x):
+            return np.exp(1j * np.angle(x)) * self.simplex(np.abs(x))
+        else:
+            return np.sign(x) * self.simplex(np.abs(x))
+

pyproximal/proximal/L0.py

@@ -3,6 +3,92 @@
 from pyproximal.ProxOperator import _check_tau
 from pyproximal.projection import L0BallProj
 from pyproximal import ProxOperator
+from pyproximal.proximal.L1 import _current_sigma
+
+
+def _hardthreshold(x, thresh):
+    r"""Hard thresholding.
+
+    Applies hard thresholding to vector ``x`` (equal to the proximity
+    operator for :math:`\|\mathbf{x}\|_0`) as shown in [1]_.
+
+    .. [1] Chen, F., Shen, L., Suter, B.W., "Computing the proximity
+       operator of the Lp norm with 0 < p < 1",
+       IET Signal Processing, 10, 2016.
+
+    Parameters
+    ----------
+    x : :obj:`numpy.ndarray`
+        Vector
+    thresh : :obj:`float`
+        Threshold
+
+    Returns
+    -------
+    x1 : :obj:`numpy.ndarray`
+        Tresholded vector
+
+    """
+    x1 = x.copy()
+    x1[np.abs(x) <= thresh] = 0
+    return x1
+
+
+class L0(ProxOperator):
+    r"""L0 norm proximal operator.
+
+    Proximal operator of the :math:`\ell_0` norm:
+    :math:`\sigma\|\mathbf{x}\|_0 = \text{count}(x_i \ne 0)`.
+
+    Parameters
+    ----------
+    sigma : :obj:`float` or :obj:`list` or :obj:`np.ndarray` or :obj:`func`, optional
+        Multiplicative coefficient of L1 norm. This can be a constant number, a list
+        of values (for multidimensional inputs, acting on the second dimension) or
+        a function that is called passing a counter which keeps track of how many
+        times the ``prox`` method has been invoked before and returns a scalar (or a list of)
+        ``sigma`` to be used.
+
+    Notes
+    -----
+    The :math:`\ell_0` proximal operator is defined as:
+
+    .. math::
+
+        \prox_{\tau \sigma \|\cdot\|_0}(\mathbf{x}) =
+        \operatorname{hard}(\mathbf{x}, \tau \sigma) =
+        \begin{cases}
+        x_i, & x_i < -\tau \sigma \\
+        0, & -\tau\sigma \leq x_i \leq \tau\sigma \\
+        x_i,  & x_i > \tau\sigma\\
+        \end{cases}
+
+    where :math:`\operatorname{hard}` is the so-called called *hard thresholding*.
+
+    """
+    def __init__(self, sigma=1.):
+        super().__init__(None, False)
+        self.sigma = sigma
+        self.count = 0
+
+    def __call__(self, x):
+        sigma = _current_sigma(self.sigma, self.count)
+        return np.sum(np.abs(x) > sigma)
+
+    def _increment_count(func):
+        """Increment counter
+        """
+        def wrapped(self, *args, **kwargs):
+            self.count += 1
+            return func(self, *args, **kwargs)
+        return wrapped
+
+    @_increment_count
+    @_check_tau
+    def prox(self, x, tau):
+        sigma = _current_sigma(self.sigma, self.count)
+        x = _hardthreshold(x, tau * sigma)
+        return x
 
 
 class L0Ball(ProxOperator):

pyproximal/proximal/L1.py

@@ -67,7 +67,7 @@ class L1(ProxOperator):
         \operatorname{soft}(\mathbf{x}, \tau \sigma) =
         \begin{cases}
         x_i + \tau \sigma, & x_i - g_i < -\tau \sigma \\
-        g_i, & -\sigma \leq x_i - g_i \leq \tau\sigma \\
+        g_i, & -\tau\sigma \leq x_i - g_i \leq \tau\sigma \\
         x_i - \tau\sigma,  & x_i - g_i > \tau\sigma\\
         \end{cases}
 

pyproximal/proximal/__init__.py

@@ -10,6 +10,7 @@
     AffineSet	                    Affines set indicator
     Quadratic                       Quadratic function
     Nonlinear	                    Nonlinear function
+    L0                              L0 Norm
     L0Ball                          L0 Ball
     L1	                            L1 Norm
     L1Ball                          L1 Ball
@@ -57,7 +58,7 @@
 from .SingularValuePenalty import *
 
 __all__ = ['Box', 'Simplex', 'Intersection', 'AffineSet', 'Quadratic',
-           'Euclidean', 'EuclideanBall', 'L0Ball', 'L1', 'L1Ball', 'L2',
+           'Euclidean', 'EuclideanBall', 'L0', 'L0Ball', 'L1', 'L1Ball', 'L2',
            'L2Convolve', 'L21', 'L21_plus_L1', 'Huber', 'Nuclear',
            'NuclearBall', 'Orthogonal', 'VStack', 'Nonlinear', 'SCAD',
            'Log', 'ETP', 'Geman', 'QuadraticEnvelopeCard', 'SingularValuePenalty']