Implement quadratic envelope of cardinality funtion

Author: marcusvaltonen

docs/source/api/index.rst

@@ -79,6 +79,7 @@ Operators
     Log
     ETP
     Geman
+    QuadraticEnvelopeCard
 
 
 Other operators

examples/plot_concave_penalties.py

@@ -55,3 +55,9 @@ def compare_penalty_and_proximal_operator(penalty):
 # The Geman penalty
 geman = pyproximal.Geman(3, 1.2)
 compare_penalty_and_proximal_operator(geman)
+
+
+###############################################################################
+# The quadratic envelope of the l0-penalty
+f_mu = pyproximal.QuadraticEnvelopeCard(1.5)
+compare_penalty_and_proximal_operator(f_mu)

pyproximal/proximal/QuadraticEnvelope.py

@@ -0,0 +1,88 @@
+import numpy as np
+
+from pyproximal.ProxOperator import _check_tau
+from pyproximal import ProxOperator
+
+
+class QuadraticEnvelopeCard(ProxOperator):
+    r"""Quadratic envelope of the :math:`\ell_0`-penalty.
+
+    The :math:`\ell_0`-penalty is also known as the *cardinality function*, and the
+    quadratic envelope :math:`\mathcal{Q}(\mu\|\cdot\|_0)` of it is defined as
+
+    .. math::
+
+        \mathcal{Q}(\mu\|\cdot\|_0)(x) = \sum_i \left(\mu - \frac{1}{2}\max(0, \sqrt{2\mu} - |x_i|)^2\right)
+
+    where :math:`\mu \geq 0`.
+
+    Parameters
+    ----------
+    mu : :obj:`float`
+        Threshold parameter.
+
+    Notes
+    -----
+    The terminology *quadratic envelope* was coined in [1]_, however, the rationale has
+    been used earlier, e.g. in [2]_. In a general setting, the quadratic envelope
+    :math:`\mathcal{Q}(f)(x)` is defined such that
+
+    .. math::
+
+        \left(f(x) + \frac{1}{2}\|x-y\|_2^2\right)^{**} = \mathcal{Q}(f)(x) + \frac{1}{2}\|x-y\|_2^2
+
+    where :math:`g^{**}` denotes the bi-conjugate of :math:`g`, which is the l.s.c.
+    convex envelope of :math:`g`.
+
+    There is no closed-form expression for :math:`\mathcal{Q}(f)(x)` given an arbitrary
+    function :math:`f`. However, for certain special cases, such as in the case of the
+    cardinality function, such expressions do exist.
+
+    The proximal operator is given by
+
+    .. math::
+
+        \prox_{\tau\mathcal{Q}(\mu\|\cdot\|_0)}(x) =
+        \begin{cases}
+        x_i, & |x_i| \geq \sqrt{2 \mu} \\
+        \frac{x_i-\tau\sqrt{2\mu}\sgn(x_i)}{1-\tau}, & \tau\sqrt{2\mu} < |x_i| < \sqrt{2 \mu} \\
+        0, & |x_i| \leq \tau\sqrt{2 \mu}
+        \end{cases}
+
+    By inspecting the structure of the proximal operator it is clear that large values
+    are unaffected, whereas smaller ones are penalized partially or completely. Such
+    properties are desirable to counter the effect of *shrinking bias* observed with
+    e.g. the :math:`\ell_1`-penalty. Note that in the limit :math:`\tau=1` this becomes
+    the hard thresholding with threshold :math:`\sqrt{2\mu}`. It should also be noted
+    that this proximal operator is identical to the Minimax Concave Penalty (MCP)
+    proposed in [3]_.
+
+    .. [1] Carlsson, M. "On Convex Envelopes and Regularization of Non-convex
+        Functionals Without Moving Global Minima", In Journal of Optimization Theory
+        and Applications, 183:66–84, 2019.
+    .. [2] Larsson, V. and Olsson, C. "Convex Low Rank Approximation", In International
+        Journal of Computer Vision (IJCV), 120:194–214, 2016.
+    .. [3] Zhang et al. "Nearly unbiased variable selection under minimax concave
+        penalty", In the Annals of Statistics, 38(2):894–942, 2010.
+
+    """
+
+    def __init__(self, mu):
+        super().__init__(None, False)
+        self.mu = mu
+
+    def __call__(self, x):
+        return np.sum(self.elementwise(x))
+
+    def elementwise(self, x):
+        return self.mu - 0.5 * np.maximum(0, np.sqrt(2 * self.mu) - np.abs(x)) ** 2
+
+    @_check_tau
+    def prox(self, x, tau):
+        r = np.abs(x)
+        idx = r < np.sqrt(2 * self.mu)
+        if tau >= 1:
+            r[idx] = 0
+        else:
+            r[idx] = np.maximum(0, (r[idx] - tau * np.sqrt(2 * self.mu)) / (1 - tau))
+        return r * np.sign(x)

pyproximal/proximal/__init__.py

@@ -28,6 +28,7 @@
     Log                             Logarithmic
     ETP                             Exponential-type penalty
     Geman                           Geman penalty
+    QuadraticEnvelopeCard           The quadratic envelope of the cardinality function
 
 """
 
@@ -51,9 +52,10 @@
 from .Log import *
 from .ETP import *
 from .Geman import *
+from .QuadraticEnvelope import *
 
 __all__ = ['Box', 'Simplex', 'Intersection', 'AffineSet', 'Quadratic',
            'Euclidean', 'EuclideanBall', 'L0Ball', 'L1', 'L1Ball', 'L2',
            'L2Convolve', 'L21', 'L21_plus_L1', 'Huber', 'Nuclear',
            'NuclearBall', 'Orthogonal', 'VStack', 'Nonlinear', 'SCAD',
-           'Log', 'ETP', 'Geman']
+           'Log', 'ETP', 'Geman', 'QuadraticEnvelopeCard']

pytests/test_concave_penalties.py

@@ -3,10 +3,10 @@
 import numpy as np
 
 from pyproximal.utils import moreau
-from pyproximal.proximal import ETP, Geman, Log, SCAD
+from pyproximal.proximal import ETP, Geman, Log, SCAD, QuadraticEnvelopeCard
 
-par1 = {'nx': 10, 'sigma': 1., 'a': 2.1, 'gamma': 0.5, 'dtype': 'float32'}  # even float32
-par2 = {'nx': 11, 'sigma': 2., 'a': 3.7, 'gamma': 5.0, 'dtype': 'float64'}  # odd float64
+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
 
 
 @pytest.mark.parametrize("par", [(par1), (par2)])
@@ -79,3 +79,19 @@ def test_Geman(par):
     # Check proximal operator
     tau = 2.
     assert moreau(geman, x, tau)
+
+
+@pytest.mark.parametrize("par", [(par1), (par2)])
+def test_QuadraticEnvelopeCard(par):
+    """QuadraticEnvelopeCard penalty and proximal/dual proximal
+    """
+    np.random.seed(10)
+    fmu = QuadraticEnvelopeCard(mu=par['mu'])
+    # Quadratic envelope of the l0-penalty
+    x = np.random.normal(0., 10.0, par['nx']).astype(par['dtype'])
+    expected = np.linalg.norm(par['mu'] - 0.5 * np.maximum(0, np.sqrt(2 * par['mu']) - np.abs(x)) ** 2, 1)
+    assert fmu(x) == pytest.approx(expected)
+
+    # Check proximal operator
+    tau = 0.25
+    assert moreau(fmu, x, tau)