Merge pull request #191 from mrava87/feat-hubercirc

feat: added HuberCircular

Author: mrava87

docs/source/api/index.rst

@@ -66,6 +66,7 @@ Convex
     EuclideanBall
     Hankel
     Huber
+    HuberCircular
     Intersection
     L0
     L0Ball

pyproximal/optimization/primal.py

@@ -345,7 +345,7 @@ def AcceleratedProximalGradient(proxf, proxg, x0, tau=None, beta=0.5,
     r"""Accelerated Proximal gradient
 
     This is a thin wrapper around :func:`pyproximal.optimization.primal.ProximalGradient` with
-    ``vandenberghe`` or ``fista``acceleration. See :func:`pyproximal.optimization.primal.ProximalGradient`
+    ``vandenberghe`` or ``fista`` acceleration. See :func:`pyproximal.optimization.primal.ProximalGradient`
     for details.
 
     """

pyproximal/proximal/Huber.py

@@ -1,14 +1,78 @@
 import numpy as np
 
 from pyproximal.ProxOperator import _check_tau
-from pyproximal.projection import BoxProj
 from pyproximal import ProxOperator
+from pyproximal.proximal import L2, L1
 
 
 class Huber(ProxOperator):
     r"""Huber norm proximal operator.
 
-    Proximal operator of the Huber norm defined as:
+    Proximal operator of the Huber norm defined as 
+    :math:`H_\alpha(\mathbf{x}) = \sum_i H_\alpha(x_i)` where:
+
+    .. math::
+
+        H_\alpha(x_i) = 
+        \begin{cases}
+        \frac{|x_i|^2}{2 \alpha}, & |x_i| \leq \alpha \\
+        |x_i| - \frac{\alpha}{2}, & |x_i| > \alpha
+        \end{cases}
+
+    which behaves like a :math:`\ell_2` norm for :math:`|x_i| \leq \alpha` and a
+    :math:`\ell_1` norm for :math:`|x_i| > \alpha`.
+
+    Parameters
+    ----------
+    alpha : :obj:`float`
+        Huber parameter
+
+    Notes
+    -----
+    The Huber proximal operator is defined as:
+
+    .. math::
+
+        \prox_{\tau H_\alpha(\cdot)}(\mathbf{x}) =
+        \begin{cases}
+        \prox_{\frac{\tau}{2 \alpha} |x_i|^2}(x_i), & |x_i| \leq \alpha \\
+        \prox_{\tau |x_i|}(x_i), & |x_i| > \alpha
+        \end{cases}
+        
+    """
+    def __init__(self, alpha):
+        super().__init__(None, False)
+        self.alpha = alpha
+        self.l2 = L2(sigma=1. / self.alpha)
+        self.l1 = L1()
+
+    def __call__(self, x):
+        h = np.zeros_like(x)
+        xabs = np.abs(x)
+        mask = xabs > self.alpha
+        h[~mask] = xabs[~mask] ** 2 / (2. * self.alpha)
+        h[mask] = xabs[mask] - self.alpha / 2.
+        return np.sum(h)
+
+    @_check_tau
+    def prox(self, x, tau):
+        y = np.zeros_like(x)
+        xabs = np.abs(x)
+        mask = xabs > self.alpha
+        y[~mask] = self.l2.prox(x[~mask], tau)
+        y[mask] = self.l1.prox(x[mask], tau)
+        # alternative from https://math.stackexchange.com/questions/1650411/
+        # proximal-operator-of-the-huber-loss-function... currently commented
+        # as it does not provide the same result
+        # y = (1. - tau / np.maximum(np.abs(x), tau + self.alpha)) * x
+
+        return y
+    
+
+class HuberCircular(ProxOperator):
+    r"""Circular Huber norm proximal operator.
+
+    Proximal operator of the Circular Huber norm defined as:
 
     .. math::
 
@@ -18,8 +82,8 @@ class Huber(ProxOperator):
         \|\mathbf{x}\|_2 - \frac{\alpha}{2}, & \|\mathbf{x}\|_2 > \alpha \\
         \end{cases}
 
-    which behaves like a :math:`\ell_2` norm for :math:`|x_i| \leq \alpha` and a
-    :math:`\ell_1` norm for :math:`|x_i| < \alpha`.
+    which behaves like a :math:`\ell_2` norm for :math:`\|\mathbf{x}\|_2 \leq \alpha` and a
+    :math:`\ell_1` norm for :math:`\|\mathbf{x}\|_2 > \alpha`.
 
     Parameters
     ----------
@@ -28,13 +92,16 @@ class Huber(ProxOperator):
 
     Notes
     -----
-    The Huber proximal operator is defined as:
+    The Circular Huber proximal operator is defined as [1]_:
 
     .. math::
 
         \prox_{\tau H_\alpha(\cdot)}(\mathbf{x}) =
         \left( 1 - \frac{\tau}{\max\{\|\mathbf{x}\|_2, \tau + \alpha \} } \right) \mathbf{x}
 
+    .. [1] O’Donoghue, B. and Stathopoulos, G. and Boyd, S. "A Splitting Method for Optimal Control", 
+        In the IEEE Transactions on Control Systems Technology, 2013.
+        
     """
     def __init__(self, alpha):
         super().__init__(None, False)
@@ -51,4 +118,4 @@ def __call__(self, x):
     @_check_tau
     def prox(self, x, tau):
         x = (1. - tau / max(np.linalg.norm(x), tau + self.alpha)) * x
-        return x
+        return x

pyproximal/proximal/__init__.py

@@ -22,6 +22,7 @@
     L21	                            L2,1 Norm
     L21_plus_L1	                    L2,1 + L1 mixed-norm
     Huber	                        Huber Norm
+    HuberCircular	                Circular Huber Norm
     TV                              Total Variation Norm
     RelaxedMumfordShah              Relaxed Mumford Shah Norm
     Nuclear                         Nuclear Norm
@@ -69,7 +70,7 @@
 
 __all__ = ['Box', 'Simplex', 'Intersection', 'AffineSet', 'Quadratic',
            'Euclidean', 'EuclideanBall', 'L0', 'L0Ball', 'L01Ball', 'L1', 'L1Ball', 'L2',
-           'L2Convolve', 'L21', 'L21_plus_L1', 'Huber', 'TV', 'RelaxedMumfordShah',
+           'L2Convolve', 'L21', 'L21_plus_L1', 'Huber', 'HuberCircular', 'TV', 'RelaxedMumfordShah',
            'Nuclear', 'NuclearBall', 'Orthogonal', 'VStack', 'Nonlinear', 'SCAD',
            'Log', 'Log1', 'ETP', 'Geman', 'QuadraticEnvelopeCard', 'SingularValuePenalty',
            'QuadraticEnvelopeCardIndicator', 'QuadraticEnvelopeRankL2',

pytests/test_norms.py

@@ -5,8 +5,8 @@
 
 from pylops.basicoperators import Identity, Diagonal, MatrixMult, FirstDerivative
 from pyproximal.utils import moreau
-from pyproximal.proximal import Box, Euclidean, L2, L1, L21, L21_plus_L1, \
-    Huber, Nuclear, RelaxedMumfordShah, TV
+from pyproximal.proximal import Euclidean, L2, L1, L21, L21_plus_L1, \
+    Huber, HuberCircular, Nuclear, RelaxedMumfordShah, TV
 
 par1 = {'nx': 10, 'sigma': 1., 'dtype': 'float32'}  # even float32
 par2 = {'nx': 11, 'sigma': 2., 'dtype': 'float64'}  # odd float64
@@ -182,8 +182,24 @@ def test_Huber(par):
     """
     hub = Huber(alpha=par['sigma'])
 
+    # norm
+    x = np.random.uniform(0., 0.1 * par['sigma'], par['nx']).astype(par['dtype'])
+    assert hub(x) == np.sum(x ** 2) / (2 * par['sigma'])
+
+    # prox / dualprox
+    tau = 2.
+    assert moreau(hub, x, tau)
+
+
+@pytest.mark.parametrize("par", [(par1), (par2)])
+def test_HuberCircular(par):
+    """Circular Huber norm and proximal/dual proximal
+    """
+    hub = HuberCircular(alpha=par['sigma'])
+
     # norm
     x = np.random.uniform(0., 0.1, par['nx']).astype(par['dtype'])
+    x = (0.1 * par['sigma']) * x / np.linalg.norm(x)  # to ensure that is smaller than sigma
     assert hub(x) == np.linalg.norm(x) ** 2 / (2 * par['sigma'])
 
     # prox / dualprox