Merge pull request #101 from mrava87/main

Added Log1 proximal

Author: mrava87

docs/source/api/index.rst

@@ -89,6 +89,7 @@ Non-Convex
     ETP
     Geman
     Log
+    Log1
     QuadraticEnvelopeCard
     QuadraticEnvelopeCardIndicator
     SCAD

pyproximal/proximal/Log.py

@@ -96,6 +96,76 @@ def prox(self, x, tau):
         r = np.vstack((np.zeros_like(c), (b[idx] + c) / (2 * self.gamma)))
         val = tau * self.elementwise(r) + (r - np.abs(x[idx])) ** 2 / 2
         idx_minima = np.argmin(val, axis=0)
-        out[idx] = r[idx_minima, range(r.shape[1])]
+        out[idx] = r[idx_minima, list(range(r.shape[1]))]
         out *= np.sign(x)
         return out
+
+
+class Log1(ProxOperator):
+    r"""Logarithmic penalty 2.
+
+    The logarithmic penalty (Log) is defined as
+
+    .. math::
+
+        \mathrm{Log}_{\sigma,\delta}(\mathbf{x}) = \sigma \sum_i \log(|x_i| + \delta)
+
+    where :math:`{\sigma>0}`, :math:`{\gamma>0}`.
+
+    Parameters
+    ----------
+    sigma : :obj:`float`
+        Multiplicative coefficient of Log norm.
+    delta : :obj:`float`, optional
+        Regularization parameter. Default is 1e-10.
+
+    Notes
+    -----
+    The logarithmic penalty gives rise to a log-thresholding that is
+    a smooth alternative falling in between the hard and soft thresholding.
+
+    The proximal operator is defined as [1]_:
+
+    .. math::
+
+        \prox_{\tau \sigma log}(\mathbf{x}) =
+        \begin{cases}
+        0.5 (x_i + \delta - \sqrt{(x_i-\delta)^2-2\tau \sigma}), & x_i < -x_0 \\
+        0, & -x_0 \leq x_i \leq  x_0 \\
+        0.5 (x_i - \delta + \sqrt{(x_i+\delta)^2-2\tau \sigma}), & x_i  > x_0\\
+        \end{cases}
+
+    where :math:`x_0=\sqrt{2 \tau \sigma} - \delta`.
+
+    .. [1] Malioutov, D., and Aravkin, A. "Iterative log thresholding",
+        Arxiv, 2013.
+
+    """
+
+    def __init__(self, sigma, delta=1e-10):
+        super().__init__(None, False)
+        if delta < 0:
+            raise ValueError('Variable "delta" must be positive.')
+        self.sigma = sigma
+        self.delta = delta
+
+    def __call__(self, x):
+        return np.sum(self.elementwise(x))
+
+    def elementwise(self, x):
+        return np.log(np.abs(x) + self.delta)
+
+    @_check_tau
+    def prox(self, x, tau):
+        tau1 = self.sigma * tau
+        thresh = np.sqrt(2*tau1) - self.delta
+        x1 = np.zeros_like(x, dtype=x.dtype)
+        if np.iscomplexobj(x):
+            x1[np.abs(x) > thresh] = 0.5 * np.exp(1j * np.angle(x[np.abs(x) > thresh])) * \
+                                     (np.abs(x[np.abs(x) > thresh]) - self.delta +
+                                      np.sqrt(np.abs(x[np.abs(x) > thresh] + self.delta) ** 2 - 2 * tau1))
+        else:
+            x1[x > thresh] = 0.5 * (x[x > thresh] - self.delta + np.sqrt(np.abs(x[x > thresh] + self.delta) ** 2 - 2 * tau1))
+            x1[x < -thresh] = 0.5 * (x[x < -thresh] + self.delta - np.sqrt(np.abs(x[x < -thresh] - self.delta) ** 2 - 2 * tau1))
+        return x1
+

pyproximal/proximal/__init__.py

@@ -28,6 +28,7 @@
     VStack	                        Stack of proximal operators
     SCAD                            Smoothly clipped absolute deviation
     Log                             Logarithmic
+    Log1                            Another form of logarithmic
     ETP                             Exponential-type penalty
     Geman                           Geman penalty
     QuadraticEnvelopeCard           The quadratic envelope of the cardinality function
@@ -67,6 +68,6 @@
            'Euclidean', 'EuclideanBall', 'L0', 'L0Ball', 'L1', 'L1Ball', 'L2',
            'L2Convolve', 'L21', 'L21_plus_L1', 'Huber', 'TV', 'Nuclear',
            'NuclearBall', 'Orthogonal', 'VStack', 'Nonlinear', 'SCAD',
-           'Log', 'ETP', 'Geman', 'QuadraticEnvelopeCard', 'SingularValuePenalty',
+           'Log', 'Log1', 'ETP', 'Geman', 'QuadraticEnvelopeCard', 'SingularValuePenalty',
            'QuadraticEnvelopeCardIndicator', 'QuadraticEnvelopeRankL2',
            'Hankel']

pytests/test_concave_penalties.py

@@ -3,7 +3,7 @@
 import numpy as np
 
 from pyproximal.utils import moreau
-from pyproximal.proximal import ETP, Geman, Log, SCAD, QuadraticEnvelopeCard
+from pyproximal.proximal import ETP, Geman, Log, Log1, SCAD, QuadraticEnvelopeCard
 
 par1 = {'nx': 10, 'sigma': 1., 'a': 2.1, 'gamma': 0.5, 'mu': 0.5, 'dtype': 'float32'}  # even float32
 par2 = {'nx': 11, 'sigma': 2., 'a': 3.7, 'gamma': 5.0, 'mu': 1.5, 'dtype': 'float64'}  # odd float64
@@ -49,6 +49,19 @@ def test_Log(par):
     assert moreau(log, x, tau)
 
 
+@pytest.mark.parametrize("par", [(par1), (par2)])
+def test_Log1(par):
+    """Log1 penalty and proximal/dual proximal
+    """
+    np.random.seed(10)
+    log = Log1(sigma=par['sigma'])
+
+    # Check proximal operator
+    tau = 2.
+    x = np.random.normal(0., 10.0, par['nx']).astype(par['dtype'])
+    assert moreau(log, x, tau)
+
+
 @pytest.mark.parametrize("par", [(par1), (par2)])
 def test_ETP(par):
     """Exponential-type penalty and proximal/dual proximal