minor: small changes to rMS and associated tutorial

Author: mrava87

docs/source/api/index.rst

@@ -92,6 +92,7 @@ Non-Convex
     Log1
     QuadraticEnvelopeCard
     QuadraticEnvelopeCardIndicator
+    RelaxedMumfordShah
     SCAD
 
 

pyproximal/proximal/rMS.py renamed to pyproximal/proximal/RelaxedMS.py

@@ -2,6 +2,7 @@
 
 from pyproximal.ProxOperator import _check_tau
 from pyproximal import ProxOperator
+from pyproximal.proximal.L1 import _current_sigma
 
 
 def _l2(x, thresh):
@@ -26,34 +27,27 @@ def _l2(x, thresh):
     return y
 
 
-def _current_sigma(sigma, count):
-    if not callable(sigma):
-        return sigma
-    else:
-        return sigma(count)
-
-
 def _current_kappa(kappa, count):
     if not callable(kappa):
         return kappa
     else:
         return kappa(count)
 
 
-class rMS(ProxOperator):
-    r"""relaxed Mumford-Shoh norm proximal operator.
+class RelaxedMumfordShah(ProxOperator):
+    r"""Relaxed Mumford-Shah norm proximal operator.
 
     Proximal operator of the relaxed Mumford-Shah norm:
-    :math:`\text{rMS}(x) = \min (\alpha\Vert x\Vert_2^2, \kappa).`.
+    :math:`\text{rMS}(x) = \min (\alpha\Vert x\Vert_2^2, \kappa)`.
 
     Parameters
     ----------
     sigma : :obj:`float` or :obj:`list` or :obj:`np.ndarray` or :obj:`func`, optional
-        Multiplicative coefficient of L2 norm that controls the smoothness. 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.
+        Multiplicative coefficient of L2 norm that controls the smoothness of the solutuon.
+        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.
     kappa : :obj:`float` or :obj:`list` or :obj:`np.ndarray` or :obj:`func`, optional
         Constant value in the rMS norm which essentially controls when the norm allows a jump. This can be a
         constant number, a list of values (for multidimensional inputs, acting on the second dimension) or
@@ -65,12 +59,12 @@ class rMS(ProxOperator):
 
     Notes
     -----
-    The :math:`\ell_1` proximal operator is defined as [1]_:
+    The :math:`rMS` proximal operator is defined as [1]_:
 
     .. math::
-        \text{prox}_{\text{rMS}}(x) =
+        \text{prox}_{\tau \text{rMS}}(x) =
         \begin{cases}
-        \frac{1}{1+2\alpha}x & \text{ if } & \vert x\vert \leq \sqrt{\frac{\kappa}{\alpha}(1 + 2\alpha)} \\
+        \frac{1}{1+2\tau\alpha}x & \text{ if } & \vert x\vert \leq \sqrt{\frac{\kappa}{\alpha}(1 + 2\tau\alpha)} \\
         \kappa & \text{ else }
         \end{cases}.
 
@@ -89,7 +83,7 @@ def __init__(self, sigma=1., kappa=1., g=None):
     def __call__(self, x):
         sigma = _current_sigma(self.sigma, self.count)
         kappa = _current_sigma(self.kappa, self.count)
-        return np.minimum(sigma * np.linalg.norm(x)**2, kappa)
+        return np.minimum(sigma * np.linalg.norm(x) ** 2, kappa)
 
     def _increment_count(func):
         """Increment counter
@@ -107,10 +101,3 @@ def prox(self, x, tau):
 
         x = np.where(np.abs(x) <= np.sqrt(kappa / sigma * (1 + 2 * tau * sigma)), _l2(x, tau * sigma), x)
         return x
-
-    @_check_tau
-    def proxdual(self, x, tau):
-        # x - tau * self.prox(x / tau, 1. / tau)
-        x = self._proxdual_moreau(x, tau)
-
-        return x

pyproximal/proximal/__init__.py

@@ -7,21 +7,22 @@
     Box                             Box indicator
     Simplex                         Simplex indicator
     Intersection                    Intersection indicator
-    AffineSet	                      Affines set indicator
+    AffineSet	                    Affines set indicator
     Quadratic                       Quadratic function
-    Nonlinear	                      Nonlinear function
+    Nonlinear	                    Nonlinear function
     L0                              L0 Norm
     L0Ball                          L0 Ball
     L1	                            L1 Norm
     L1Ball                          L1 Ball
-    Euclidean	                      Euclidean Norm
-    EuclideanBall	                  Euclidean Ball
+    Euclidean	                    Euclidean Norm
+    EuclideanBall	                Euclidean Ball
     L2	                            L2 Norm
     L2Convolve	                    L2 Norm of convolution operator
     L21	                            L2,1 Norm
     L21_plus_L1	                    L2,1 + L1 mixed-norm
-    Huber	                          Huber Norm
-    TV                              Total Variation Norm                                   
+    Huber	                        Huber Norm
+    TV                              Total Variation Norm
+    RelaxedMumfordShah              Relaxed Mumford Shah Norm
     Nuclear                         Nuclear Norm
     NuclearBall                     Nuclear Ball
     Orthogonal	                    Product between orthogonal operator and vector
@@ -53,6 +54,7 @@
 from .L21_plus_L1 import *
 from .Huber import *
 from .TV import *
+from .RelaxedMS import *
 from .Nuclear import *
 from .Orthogonal import *
 from .VStack import *
@@ -66,8 +68,8 @@
 
 __all__ = ['Box', 'Simplex', 'Intersection', 'AffineSet', 'Quadratic',
            'Euclidean', 'EuclideanBall', 'L0', 'L0Ball', 'L1', 'L1Ball', 'L2',
-           'L2Convolve', 'L21', 'L21_plus_L1', 'Huber', 'TV', 'Nuclear',
-           'NuclearBall', 'Orthogonal', 'VStack', 'Nonlinear', 'SCAD',
-           'Log', 'Log1', 'ETP', 'Geman', 'QuadraticEnvelopeCard', 'rMS', 'SingularValuePenalty',
+           'L2Convolve', 'L21', 'L21_plus_L1', 'Huber', 'TV', 'RelaxedMumfordShah',
+           'Nuclear', 'NuclearBall', 'Orthogonal', 'VStack', 'Nonlinear', 'SCAD',
+           'Log', 'Log1', 'ETP', 'Geman', 'QuadraticEnvelopeCard', 'SingularValuePenalty',
            'QuadraticEnvelopeCardIndicator', 'QuadraticEnvelopeRankL2',
            'Hankel']

pytests/test_norms.py

@@ -5,7 +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, TV
+from pyproximal.proximal import Box, Euclidean, L2, L1, L21, L21_plus_L1, \
+    Huber, Nuclear, RelaxedMumfordShah, TV
 
 par1 = {'nx': 10, 'sigma': 1., 'dtype': 'float32'}  # even float32
 par2 = {'nx': 11, 'sigma': 2., 'dtype': 'float64'}  # odd float64
@@ -202,6 +203,22 @@ def test_TV(par):
     assert_array_almost_equal(tv(x), par['sigma'] * np.sum(np.abs(dx), axis=0))
 
 
+@pytest.mark.parametrize("par", [(par1), (par2)])
+def test_rMS(par):
+    """rMS norm and proximal/dual proximal
+    """
+    kappa = 1.
+    rMS = RelaxedMumfordShah(sigma=par['sigma'], kappa=kappa)
+
+    # norm
+    x = np.random.normal(0., 1., par['nx']).astype(par['dtype'])
+    assert rMS(x) == np.minimum(par['sigma'] * np.linalg.norm(x) ** 2, kappa)
+
+    # prox / dualprox
+    tau = 2.
+    assert moreau(rMS, x, tau)
+
+
 def test_Nuclear_FOM():
     """Nuclear norm benchmark with FOM solver
     """